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 /* */
7 /* file: write_mesh_gmv.c */
8 /* */
9 /* */
10 /* description: This program converts ALBERTA meshes and dof_real_[d_]vecs */
11 /* to a GMV 4.0 format (ASCII or "iecxi?r?", one of the binary */
12 /* format versions, NOT NECESSARILY PORTABLE) */
13 /*--------------------------------------------------------------------------*/
14 /* */
15 /* authors: Daniel Koester */
16 /* Institut fuer Mathematik */
17 /* Universitaet Augsburg */
18 /* Universitaetsstr. 14 */
19 /* D-86159 Augsburg, Germany */
20 /* */
21 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
22 /* */
23 /* (c) by D. Koester (2004) */
24 /*--------------------------------------------------------------------------*/
25
26 #ifdef HAVE_CONFIG_H
27 # include "config.h"
28 #endif
29
30 #include "alberta.h"
31 #include "alberta_intern.h"
32
33 #include <ctype.h>
34
35 #define VERT_IND(dim,i,j) ((i)*N_VERTICES(dim)+(j))
36
37 static const REAL_B center_bary[DIM_LIMIT+1] = {
38 INIT_BARY_0D(0.0),
39 INIT_BARY_1D(0.5, 0.5),
40 INIT_BARY_2D(1.0/3.0, 1.0/3.0, 1.0/3.0),
41 INIT_BARY_3D(0.25,0.25,0.25,0.25)
42 };
43
44 /***************************************************************************/
45 /* gmv_open_ascii(filenam, mesh, sim_time): open a file for ASCII GMV */
46 /* output. */
47 /***************************************************************************/
48
gmv_open_ascii(const char * filenam,MESH * mesh,REAL sim_time)49 static FILE * gmv_open_ascii(const char *filenam, MESH *mesh, REAL sim_time)
50 {
51 /* FUNCNAME("gmv_open_ascii"); */
52 FILE *file = NULL;
53
54 file = fopen(filenam, "w");
55
56 if (file == NULL) {
57 return NULL;
58 }
59
60 /* Write header. */
61 fprintf(file, "gmvinput ascii\n");
62
63 if(mesh->name) {
64 fprintf(file, "comments\n");
65 fprintf(file, "Mesh '%s'\n", mesh->name);
66 fprintf(file, "endcomm\n");
67 }
68
69 if (isfinite(sim_time)) {
70 fprintf(file, "probtime %.6E\n", sim_time);
71 }
72 fprintf(file, "codename ALBERTA \n");
73 fprintf(file, "codever 2.0 \n");
74
75 return file;
76 }
77
78 /***************************************************************************/
79 /* gmv_open_bin(filenam,isize,rsize,time): open a file for binary output */
80 /***************************************************************************/
81
gmv_open_bin(const char * filenam,int isize,int rsize,REAL sim_time)82 static FILE * gmv_open_bin(const char *filenam, int isize, int rsize,
83 REAL sim_time)
84 {
85 FUNCNAME("gmv_open_bin");
86 FILE *fp = NULL;
87
88 /* Check if this machine can address 64bit integers (isize = 8). */
89 if (isize == 8)
90 TEST_EXIT (sizeof(long) >= 8,
91 "Cannot write 8 byte integers on this machine.\n");
92
93 fp = fopen(filenam, "w");
94
95 if (fp == NULL) {
96 return NULL;
97 }
98
99 /* Write header. */
100 AI_fwrite("gmvinput", sizeof(char), 8, fp);
101 if (isize == 4 && rsize == 8) {
102 AI_fwrite("iecxi4r8", sizeof(char), 8, fp);
103 }
104 else if (isize == 8 && rsize == 4) {
105 AI_fwrite("iecxi8r4", sizeof(char), 8, fp);
106 }
107 else if (isize == 8 && rsize == 8) {
108 AI_fwrite("iecxi8r8", sizeof(char), 8, fp);
109 }
110 else /* just do the normal "iecxi4r4" thing */ {
111 AI_fwrite("iecxi4r4", sizeof(char), 8, fp);
112 }
113
114 AI_fwrite("probtime", sizeof(char), 8, fp);
115 AI_fwrite(&sim_time, sizeof(REAL), 1, fp);
116 AI_fwrite("codenameALBERTA codever 2.0 ", sizeof(char), 32, fp);
117
118 return fp;
119 }
120
121
122 /***************************************************************************/
123 /* convert_string(in): Replace all white-space characters in "string" */
124 /* with underscores. */
125 /***************************************************************************/
126
convert_string(char * pos)127 static void convert_string(char *pos)
128 {
129 for(; *pos; pos++)
130 if(!isgraph((int)*pos)) *pos = '_';
131 }
132
133
134 /***************************************************************************/
135 /* write_coord_array(file, n_vert, coords, write_ascii): dump the given */
136 /* coordinate array to the GMV file. */
137 /***************************************************************************/
138
write_coord_array(FILE * file,int n_vert,REAL_D * coords,bool write_ascii)139 static void write_coord_array(FILE *file, int n_vert, REAL_D *coords,
140 bool write_ascii)
141 {
142 REAL *tempxyz;
143 int i, j;
144
145 if(write_ascii) {
146 fprintf(file, "nodev %d\n", n_vert);
147
148 for(i = 0; i < n_vert; i++) {
149 for(j = 0; j < 3; j++) {
150 if(j < DIM_OF_WORLD)
151 fprintf(file, "%.10E ", coords[i][j]);
152 else
153 fprintf(file, "0.0 ");
154 }
155 fprintf(file, "\n");
156 }
157 }
158 else {
159 AI_fwrite("nodev ",sizeof(char),8, file);
160
161 AI_fwrite(&n_vert, sizeof(int), 1, file);
162
163 #if DIM_OF_WORLD == 3
164 tempxyz = (REAL *) coords;
165 #else
166 tempxyz = MEM_CALLOC(n_vert * 3, REAL);
167
168 for (i = 0; i < n_vert; i++)
169 for(j = 0; j < DIM_OF_WORLD; j++)
170 tempxyz[3*i + j] = coords[i][j];
171 #endif
172
173 AI_fwrite(tempxyz, sizeof(REAL), 3*n_vert, file);
174
175 #if DIM_OF_WORLD < 3
176 MEM_FREE(tempxyz, n_vert * 3, REAL);
177 #endif
178 }
179
180 return;
181 }
182
183
184 /***************************************************************************/
185 /* write_elem_array(file, dim, n_elem, elem, write_ascii): dump the given */
186 /* element array to the GMV file. */
187 /***************************************************************************/
188
write_elem_array(FILE * file,int dim,int n_elem,int * elem,bool write_ascii)189 static void write_elem_array(FILE *file, int dim, int n_elem, int *elem,
190 bool write_ascii)
191 {
192 const char *cell_types[4] = {"general ", "line ", "tri ", "tet "};
193 int i, j, dim_1 = dim + 1, one = 1, tmp;
194
195 if(write_ascii) {
196 fprintf(file, "cells %d\n", n_elem);
197
198 for(i = 0; i < n_elem; i++) {
199 fprintf(file, "%s%d\n", cell_types[dim], dim+1);
200
201 for(j = 0; j < N_VERTICES(dim); j++) {
202 if(dim == 0)
203 fprintf(file, "1 ");
204
205 fprintf(file, "%d ", 1 + elem[VERT_IND(dim,i,j)]);
206 }
207 fprintf(file, "\n");
208 }
209 }
210 else {
211 AI_fwrite("cells ", sizeof(char), 8, file);
212 AI_fwrite(&n_elem, sizeof(int), 1, file);
213
214 for(i = 0; i < n_elem; i++) {
215 AI_fwrite(cell_types[dim], sizeof(char), 8, file);
216 AI_fwrite(&dim_1, sizeof(int), 1, file);
217
218 for(j = 0; j < N_VERTICES(dim); j++) {
219 if(dim == 0)
220 AI_fwrite(&one, sizeof(int), 1, file);
221
222 tmp = 1 + elem[VERT_IND(dim, i, j)];
223 AI_fwrite(&tmp, sizeof(int), 1, file);
224 }
225 }
226 }
227
228 return;
229 }
230
231
232 /***************************************************************************/
233 /* write_mat_array(file, n_elem, mat_array, write_ascii): dump the given */
234 /* material array to the GMV file. */
235 /***************************************************************************/
236
write_mat_array(FILE * file,int n_elem,int * mat,bool write_ascii)237 static void write_mat_array(FILE *file, int n_elem, int *mat, bool write_ascii)
238 {
239 const char *mat_types[2] = {"affine ", "parametr"};
240 int i, tmp;
241
242 if(write_ascii) {
243 fprintf(file, "material 2 0\n");
244 fprintf(file, "%s\n", mat_types[0]);
245 fprintf(file, "%s\n", mat_types[1]);
246
247 for(i = 0; i < n_elem; i++)
248 fprintf(file, "%d\n", mat[i]);
249 }
250 else {
251 AI_fwrite("material", sizeof(char), 8, file);
252 tmp = 2;
253 AI_fwrite(&tmp, sizeof(int), 1, file);
254 tmp = 0;
255 AI_fwrite(&tmp, sizeof(int), 1, file);
256 AI_fwrite(mat_types[0], sizeof(char), 32, file);
257 AI_fwrite(mat_types[1], sizeof(char), 32, file);
258
259 AI_fwrite(mat, sizeof(int), n_elem, file);
260 }
261
262 return;
263 }
264
265
266 /***************************************************************************/
267 /* add_refined_data(file, mesh, use_refined_mesh, write_ascii, */
268 /* write_mesh_date, n_drv, */
269 /* drv_ptr, n_drv_d, drv_d_ptr, velocity): */
270 /* */
271 /* This routine enables us to output more information for higher order */
272 /* Lagrange finite elements. Since GMV only displays linear data we first */
273 /* perform a virtual refinement of all elements, see the array */
274 /* "n_sub_elements" below. The refined mesh is then output to GMV as usual.*/
275 /***************************************************************************/
276
277 static const int n_sub_elements[/*dimension*/ 4][/*degree*/ 5] = {
278 {1,1,1,1,1},
279 {1,1,2,3,4},
280 {1,1,4,9,16},
281 {1,1,4,10,20}
282 };
283
284 /* sub_elements[dimension][degree][max_n_elements == 20][vertex] */
285 static const int sub_elements[4][5][20][4] =
286 {/*dim=0*/ {/*p=0*/ {{0}},
287 /*p=1*/ {{0}},
288 /*p=2*/ {{0}},
289 /*p=3*/ {{0}},
290 /*p=4*/ {{0}}},
291 /*dim=1*/ {/*p=0*/ {{0,1}},
292 /*p=1*/ {{0,1}},
293 /*p=2*/ {{0,2}, {2,1}},
294 /*p=3*/ {{0,2}, {2,3}, {3,1}},
295 /*p=4*/ {{0,2}, {2,3}, {3,4}, {4,1}}},
296 /*dim=2*/ {/*p=0*/ {{0,1,2}},
297 /*p=1*/ {{0,1,2}},
298 /*p=2*/ {{0,5,4},{5,3,4},
299 {5,1,3},{4,3,2}},
300 /*p=3*/ {{0,7,6},{7,9,6},{7,8,9},
301 {6,9,5},{8,3,9},{9,4,5},
302 {8,1,3},{9,3,4},{5,4,2}},
303 /*p=4*/ {{ 0, 9, 8},{ 9,12, 8},{ 9,10,12},{ 8,12, 7},
304 {10,13,12},{12,14, 7},{10,11,13},{12,13,14},
305 { 7,14, 6},{11, 3,13},{13, 4,14},{14, 5, 6},
306 {11, 1, 3},{13, 3, 4},{14, 4, 5},{ 6, 5, 2}}},
307 /*dim=3*/ {/*p=0*/ {{ 0, 1, 2, 3}},
308 /*p=1*/ {{ 0, 1, 2, 3}},
309 /*p=2*/ {{ 0, 4, 5, 6},{ 4, 1, 7, 8},
310 { 5, 7, 2, 9},
311 { 6, 8, 9, 3}},
312 /*p=3*/ {{ 0, 4, 6, 8},{ 4, 5,19,18},{ 5, 1,10,12},
313 { 6,19, 7,17},{19,10,11,16},
314 { 7,11, 2,14},
315 { 8,18,17, 9},{18,12,16,13},
316 {17,16,14,15},
317 { 9,13,15, 3}},
318 /*p=4*/ {{ 0, 4, 7,10},{ 4, 5,31,28},{ 5, 6,32,29},{ 6, 1,13,16},
319 { 7,31, 8,25},{31,32,33,34},{32,13,14,22},
320 { 8,33, 9,26},{33,14,15,23},
321 { 9,15, 2,19},
322 {10,28,25,11},{28,29,34,30},{29,16,22,17},
323 {25,34,26,27},{34,22,23,24},
324 {26,23,19,20},
325 {11,30,27,12},{30,17,24,18},
326 {27,24,20,21},
327 {12,18,21, 3}}}};
328
identity_loc(REAL_D world,const EL_INFO * el_info,const QUAD * quad,int iq,void * dummy)329 static const REAL *identity_loc(REAL_D world, const EL_INFO *el_info,
330 const QUAD *quad, int iq, void *dummy)
331 {
332 static REAL_D space;
333
334 if (world == NULL) {
335 world = space;
336 }
337
338 if ((el_info->fill_flag & FILL_COORDS) == 0) {
339 const PARAMETRIC *parametric = el_info->mesh->parametric;
340 REAL_D world;
341
342 parametric->coord_to_world(el_info, NULL, 1, quad->lambda + iq,
343 (REAL_D *)world);
344 } else {
345 coord_to_world(el_info, quad->lambda[iq], world);
346 }
347
348 return world;
349 }
350
add_refined_data(FILE * file,MESH * mesh,bool use_refined_mesh,bool write_ascii,bool write_mesh_data,const int n_drv,DOF_REAL_VEC ** drv_ptr,const int n_drv_d,DOF_REAL_VEC_D ** drv_d_ptr,DOF_REAL_VEC_D * velocity)351 static int add_refined_data(FILE *file, MESH *mesh,
352 bool use_refined_mesh,
353 bool write_ascii,
354 bool write_mesh_data,
355 const int n_drv, DOF_REAL_VEC **drv_ptr,
356 const int n_drv_d, DOF_REAL_VEC_D **drv_d_ptr,
357 DOF_REAL_VEC_D *velocity)
358 {
359 FUNCNAME("add_refined_data");
360 int count, max_degree = -1, ne, nv, n_bas_fcts, i, j, k, n_sub;
361 int dim = mesh->dim, index, node_m = 0, n0;
362 DOF_REAL_VEC *drv = NULL;
363 DOF_REAL_VEC_D *drv_d = NULL;
364 const FE_SPACE *max_fe_space = NULL;
365 const BAS_FCTS *max_bas_fcts;
366 const BAS_FCTS *lag, *bas_fcts;
367 DOF_INT_VEC *dof_vert_ind;
368 DOF *local_dofs = NULL;
369 DOF **local_dof_ptr;
370 int *vert_ind = NULL, *elements;
371 REAL_D *vertices, buffer;
372 int *mat_array = NULL;
373 REAL *new_vecs[250] = { NULL, };
374 REAL *new_vecs_d[250][DIM_OF_WORLD] = { { NULL, }, };
375 REAL *new_vel[3] = { 0, };
376 REAL *local_new_vec;
377 REAL_D *local_new_vec_d;
378 const BAS_FCTS *new_bas_fcts;
379 const DOF_ADMIN *new_admin;
380 char name[32] = { '\0', };
381 bool kill_fe_space = false;
382 const REAL_B *lag_nodes;
383
384 /***************************************************************************/
385 /* Step 1: Look for the largest Lagrange degree of vectors on mesh. */
386 /* Remember the corresponding FE_SPACE, since we need a mapping from local */
387 /* DOFs to global indices. */
388 /***************************************************************************/
389
390 for(count = 0; count < n_drv; count++) {
391
392 drv = drv_ptr[count];
393
394 if(!drv) {
395 ERROR_EXIT("Could not find DOF_REAL_VEC %d!\n", count);
396 }
397 bas_fcts = drv->fe_space->bas_fcts;
398 lag = get_lagrange(bas_fcts->dim, bas_fcts->degree);
399 if (lag != bas_fcts) {
400 if(use_refined_mesh) {
401 WARNING("Refinement only implemented for Lagrange Finite Elements!\n");
402 WARNING("Using standard output mechanism.\n");
403 }
404 max_fe_space = drv->fe_space;
405 max_degree = 0;
406 break;
407 }
408 if(drv->fe_space->bas_fcts->degree > max_degree) {
409 max_fe_space = drv->fe_space;
410 max_degree = max_fe_space->bas_fcts->degree;
411 }
412 }
413 for(count = 0; count < n_drv_d; count++) {
414 drv_d = drv_d_ptr[count];
415
416 if(!drv_d) {
417 ERROR_EXIT("Could not find DOF_REAL_D_VEC %d!\n", count);
418 }
419 bas_fcts = drv_d->fe_space->bas_fcts;
420 lag = get_lagrange(bas_fcts->dim, bas_fcts->degree);
421 if (lag != bas_fcts) {
422 if(use_refined_mesh) {
423 WARNING("Refinement only implemented for Lagrange Finite Elements!\n");
424 WARNING("Using standard output mechanism.\n");
425 }
426 max_fe_space = drv_d->fe_space;
427 max_degree = 0;
428 break;
429 }
430 if(drv_d->fe_space->bas_fcts->degree > max_degree) {
431 max_fe_space = drv_d->fe_space;
432 max_degree = max_fe_space->bas_fcts->degree;
433 }
434 }
435 if(velocity) {
436 bas_fcts = velocity->fe_space->bas_fcts;
437 lag = get_lagrange(bas_fcts->dim, bas_fcts->degree);
438 if (lag != bas_fcts) {
439 if(use_refined_mesh) {
440 WARNING("Only implemented for Lagrange Finite Elements!\n");
441 WARNING("Using standard output mechanism.\n");
442 }
443 max_fe_space = velocity->fe_space;
444 max_degree = 0;
445 }
446
447 if(velocity->fe_space->bas_fcts->degree > max_degree) {
448 max_fe_space = velocity->fe_space;
449 max_degree = max_fe_space->bas_fcts->degree;
450 }
451 }
452
453 if(max_degree > 1) {
454 if(use_refined_mesh)
455 MSG("Maximal degree found:%d\n", max_degree);
456 else
457 max_degree = 1;
458 }
459
460 /* No vectors at all? Or no vectors with VERTEX DOFs? Then define */
461 /* max_fe_space as a new FE_SPACE having VERTEX DOFs as well as standard */
462 /* linear Lagrange basis functions. */
463
464 if(!max_fe_space || max_degree < 1) {
465 const BAS_FCTS *lagrange = get_lagrange(mesh->dim, 1);
466
467 max_fe_space = get_fe_space(mesh, "vertex FE_SPACE",
468 lagrange, -1, ADM_FLAGS_DFLT);
469 kill_fe_space = true;
470 max_degree = 1;
471 }
472
473
474 /***************************************************************************/
475 /* Step 2: Count the number of needed vertices and elements. Fill element */
476 /* and vertex arrays. */
477 /***************************************************************************/
478 n_sub = n_sub_elements[dim][max_degree];
479 n_bas_fcts = max_fe_space->bas_fcts->n_bas_fcts;
480 local_dofs = MEM_ALLOC(n_bas_fcts, DOF);
481 max_bas_fcts = max_fe_space->bas_fcts;
482 dof_vert_ind = get_dof_int_vec("vertex indices", max_fe_space);
483
484 GET_DOF_VEC(vert_ind, dof_vert_ind);
485 FOR_ALL_DOFS(max_fe_space->admin, vert_ind[dof] = -1);
486
487 elements = MEM_ALLOC(mesh->n_elements * n_sub * N_VERTICES(dim), int);
488
489 ne = nv = 0;
490 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL) {
491
492 GET_DOF_INDICES(max_bas_fcts, el_info->el, max_fe_space->admin, local_dofs);
493
494 for(i = 0; i < n_bas_fcts; i++)
495 if (vert_ind[local_dofs[i]] == -1) {
496 /* Assign a global index to each vertex. */
497 vert_ind[local_dofs[i]] = nv;
498
499 nv++;
500 }
501
502 /* Assign element indices. */
503 for(i = 0; i < n_sub; i++) {
504 for(j = 0; j < N_VERTICES(dim); j++)
505 elements[VERT_IND(dim,ne,j)] =
506 vert_ind[local_dofs[sub_elements[dim][max_degree][i][j]]];
507
508 ne++;
509 }
510 } TRAVERSE_NEXT();
511
512 /* Do a quick check on the element number. */
513 TEST_EXIT(ne == mesh->n_elements * n_sub,
514 "Wrong no. of elements counted!\n");
515
516 /***************************************************************************/
517 /* Step 3: Allocate and fill the vertex array, as well as the DOF_REAL_VEC */
518 /* and DOF_REAL_D_VEC arrays. */
519 /***************************************************************************/
520 vertices = MEM_ALLOC(nv, REAL_D);
521 local_new_vec = MEM_ALLOC(n_bas_fcts, REAL);
522 local_new_vec_d = MEM_ALLOC(n_bas_fcts, REAL_D);
523
524 for(i = 0; i < n_drv; i++) {
525 if(!drv_ptr[i]) {
526 ERROR("Could not find DOF_REAL_VEC %d!\n", i);
527 break;
528 }
529 if(drv_ptr[i]->fe_space->bas_fcts->degree == 0)
530 new_vecs[i] = MEM_ALLOC(ne, REAL);
531 else
532 new_vecs[i] = MEM_ALLOC(nv, REAL);
533 }
534 for(i = 0; i < n_drv_d; i++) {
535 if(!drv_d_ptr[i]) {
536 ERROR("Could not find DOF_REAL_D_VEC %d!\n", i);
537 break;
538 }
539 if(drv_d_ptr[i]->fe_space->bas_fcts->degree == 0)
540 index = ne;
541 else
542 index = nv;
543
544 for(j = 0; j < DIM_OF_WORLD; j++)
545 new_vecs_d[i][j] = MEM_ALLOC(index, REAL);
546 }
547 if(velocity) {
548 if(velocity->fe_space->bas_fcts->degree == 0)
549 index = ne;
550 else
551 index = nv;
552 for(i = 0; i < 3; i++)
553 new_vel[i] = MEM_CALLOC(index, REAL);
554 }
555
556 lag_nodes = LAGRANGE_NODES(max_fe_space->bas_fcts);
557
558 index = 0;
559 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_COORDS) {
560 DEF_EL_VEC_CONST(REAL_D, local_coords, N_VERTICES_MAX, N_VERTICES_MAX);
561 index++;
562
563 GET_DOF_INDICES(max_bas_fcts, el_info->el, max_fe_space->admin, local_dofs);
564
565 if (mesh->parametric) {
566 mesh->parametric->init_element(el_info, mesh->parametric);
567 }
568
569 INTERPOL_D(max_bas_fcts,
570 local_coords, el_info, -1, 0, NULL, identity_loc, NULL);
571
572 for(i = 0; i < n_bas_fcts; i++)
573 COPY_DOW(local_coords->vec[i], vertices[vert_ind[local_dofs[i]]]);
574
575 for(i = 0; i < n_drv; i++) {
576 if(drv_ptr[i]) {
577 new_bas_fcts = drv_ptr[i]->fe_space->bas_fcts;
578 new_admin = drv_ptr[i]->fe_space->admin;
579
580 if(new_bas_fcts->degree == 0) {
581 n0 = new_admin->n0_dof[CENTER];
582 node_m = mesh->node[CENTER];
583
584 local_dof_ptr = el_info->el->dof;
585
586 for(j = 0; j < n_sub; j++)
587 new_vecs[i][index*n_sub + j] =
588 drv_ptr[i]->vec[local_dof_ptr[node_m][n0]];
589 }
590 else {
591 const EL_REAL_VEC *local_new_vec =
592 fill_el_real_vec(NULL, el_info->el, drv_ptr[i]);
593 for(j = 0; j < n_bas_fcts; j++)
594 new_vecs[i][vert_ind[local_dofs[j]]] =
595 eval_uh(lag_nodes[j], local_new_vec, new_bas_fcts);
596 }
597 }
598 }
599
600 for(i = 0; i < n_drv_d; i++) {
601 if(drv_d_ptr[i]) {
602 const EL_REAL_VEC_D *local_new_vec;
603
604 new_bas_fcts = drv_d_ptr[i]->fe_space->bas_fcts;
605 new_admin = drv_d_ptr[i]->fe_space->admin;
606 local_new_vec = fill_el_real_vec_d(NULL, el_info->el, drv_d_ptr[i]);
607
608 if(new_bas_fcts->degree == 0) {
609 eval_uh_dow(buffer, center_bary[dim], local_new_vec, new_bas_fcts);
610 for(j = 0; j < n_bas_fcts; j++) {
611 for(k = 0; k < DIM_OF_WORLD; k++) {
612 new_vecs_d[i][k][vert_ind[local_dofs[j]]] = buffer[k];
613 }
614 }
615 } else {
616 for(j = 0; j < n_bas_fcts; j++) {
617 eval_uh_dow(buffer, lag_nodes[j], local_new_vec, new_bas_fcts);
618 for(k = 0; k < DIM_OF_WORLD; k++) {
619 new_vecs_d[i][k][vert_ind[local_dofs[j]]] = buffer[k];
620 }
621 }
622 }
623 }
624 }
625
626 if(velocity) {
627 const EL_REAL_VEC_D *local_new_vec;
628
629 new_bas_fcts = velocity->fe_space->bas_fcts;
630 new_admin = velocity->fe_space->admin;
631 local_new_vec = fill_el_real_vec_d(NULL, el_info->el, velocity);
632
633 if(new_bas_fcts->degree == 0) {
634 eval_uh_dow(buffer, center_bary[dim], local_new_vec, new_bas_fcts);
635 for(i = 0; i < n_bas_fcts; i++) {
636 for(j = 0; j < DIM_OF_WORLD; j++) {
637 new_vel[j][vert_ind[local_dofs[i]]] = buffer[j];
638 }
639 }
640 } else {
641 local_new_vec = fill_el_real_vec_d(NULL, el_info->el, velocity);
642 for(i = 0; i < n_bas_fcts; i++) {
643 eval_uh_dow(buffer, lag_nodes[i], local_new_vec, new_bas_fcts);
644 for(j = 0; j < DIM_OF_WORLD; j++) {
645 new_vel[j][vert_ind[local_dofs[i]]] = buffer[j];
646 }
647 }
648 }
649 }
650
651 } TRAVERSE_NEXT();
652
653 /***************************************************************************/
654 /* Step 5: Write a material array in the parametric case. */
655 /***************************************************************************/
656
657 if(mesh->parametric) {
658 mat_array = MEM_CALLOC(ne, int);
659 ne = 0;
660
661 TRAVERSE_FIRST(mesh, -1, CALL_LEAF_EL|FILL_COORDS|FILL_PROJECTION) {
662
663 if(mesh->parametric->init_element(el_info, mesh->parametric))
664 for(i = 0; i < n_sub; i++)
665 mat_array[ne + i] = 2;
666 else
667 for(i = 0; i < n_sub; i++)
668 mat_array[ne + i] = 1;
669
670 ne += n_sub;
671 } TRAVERSE_NEXT();
672
673 ne = mesh->n_elements * n_sub;
674 }
675
676
677 /***************************************************************************/
678 /* Step 6: Dump all vectors to the file. */
679 /***************************************************************************/
680
681 if(write_mesh_data) {
682 write_coord_array(file, nv, vertices, write_ascii);
683 write_elem_array(file, dim, ne, elements, write_ascii);
684 if(mesh->parametric)
685 write_mat_array(file, ne, mat_array, write_ascii);
686 }
687
688 if(n_drv && drv_ptr) {
689 if(write_ascii)
690 fprintf(file, "variable\n");
691 else
692 AI_fwrite("variable",sizeof(char),8,file);
693 }
694
695 for(i = 0; i < n_drv; i++)
696 if(new_vecs[i]) {
697 snprintf(name, 32, "%s", drv_ptr[i]->name);
698 convert_string(name);
699
700 if(drv_ptr[i]->fe_space->bas_fcts->degree == 0) {
701 j = 0;
702 index = ne*n_sub;
703 }
704 else {
705 j = 1;
706 index = nv;
707 }
708
709 if(write_ascii) {
710 fprintf(file, "%s %d\n", name, j);
711
712 for(j = 0; j < index; j++)
713 fprintf(file, "%.10E\t", new_vecs[i][j]);
714 }
715 else {
716 AI_fwrite(name, sizeof(char), 32, file);
717 AI_fwrite(&j, sizeof(int), 1, file);
718 AI_fwrite(new_vecs[i], sizeof(REAL), index, file);
719 }
720 }
721
722 if(n_drv && drv_ptr) {
723 if(write_ascii)
724 fprintf(file, "endvars\n");
725 else
726 AI_fwrite("endvars ",sizeof(char),8,file);
727 }
728
729 if(n_drv_d && drv_d_ptr) {
730 if(write_ascii)
731 fprintf(file, "vectors \n");
732 else
733 AI_fwrite("vectors ",sizeof(char),8,file);
734 }
735
736 for(i = 0; i < n_drv_d; i++)
737 if(new_vecs_d[i][0]) {
738 snprintf(name, 32, "%s", drv_d_ptr[i]->name);
739 convert_string(name);
740
741 if(drv_d_ptr[i]->fe_space->bas_fcts->degree == 0) {
742 j = 0;
743 index = ne*n_sub;
744 }
745 else {
746 j = 1;
747 index = nv;
748 }
749
750 if(write_ascii) {
751 fprintf(file, "%s %d %d %d\n", name, j, DIM_OF_WORLD, 0);
752
753 for(j = 0; j < DIM_OF_WORLD; j++)
754 for(k = 0; k < index; k++)
755 fprintf(file, "%.10E\t", new_vecs_d[i][j][k]);
756 }
757 else {
758 AI_fwrite(name, sizeof(char), 32, file);
759 AI_fwrite(&j, sizeof(int), 1, file);
760 j = DIM_OF_WORLD;
761 AI_fwrite(&j, sizeof(int), 1, file);
762 j = 0;
763 AI_fwrite(&j, sizeof(int), 1, file);
764 for(j = 0; j < DIM_OF_WORLD; j++)
765 AI_fwrite(new_vecs_d[i][j], sizeof(REAL), index, file);
766 }
767 }
768
769 if(n_drv_d && drv_d_ptr) {
770 if(write_ascii)
771 fprintf(file, "endvect\n");
772 else
773 AI_fwrite("endvect ",sizeof(char),8,file);
774 }
775
776 if(velocity) {
777 if(velocity->fe_space->bas_fcts->degree == 0) {
778 j = 0;
779 index = ne;
780 }
781 else {
782 j = 1;
783 index = nv;
784 }
785
786 if(write_ascii) {
787 fprintf(file, "velocity %d\n", j);
788
789 for(i = 0; i < 3; i++)
790 for(j = 0; j < index; j++)
791 fprintf(file, "%.10E\t", new_vel[i][j]);
792 }
793 else {
794 AI_fwrite("velocity", sizeof(char), 8, file);
795 AI_fwrite(&j, sizeof(int), 1, file);
796
797 for(i = 0; i < 3; i++)
798 AI_fwrite(new_vel[i], sizeof(REAL), index, file);
799 }
800 }
801
802 /* Free all allocated memory. . */
803 MEM_FREE(elements, ne * N_VERTICES(dim), int);
804 MEM_FREE(vertices, nv, REAL_D);
805 MEM_FREE(local_dofs, n_bas_fcts, DOF);
806 MEM_FREE(local_new_vec, n_bas_fcts, REAL);
807 MEM_FREE(local_new_vec_d, n_bas_fcts, REAL_D);
808 if(mesh->parametric)
809 MEM_FREE(mat_array, ne, int);
810
811 for(i = 0; i < n_drv; i++)
812 if(new_vecs[i]) {
813 if(drv_ptr[i]->fe_space->bas_fcts->degree == 0)
814 MEM_FREE(new_vecs[i], ne, REAL);
815 else
816 MEM_FREE(new_vecs[i], nv, REAL);
817 }
818 for(i = 0; i < n_drv_d; i++)
819 if(new_vecs_d[i][0]) {
820 for(j = 0; j < DIM_OF_WORLD; j++)
821 if(drv_d_ptr[i]->fe_space->bas_fcts->degree == 0)
822 MEM_FREE(new_vecs_d[i][j], ne, REAL);
823 else
824 MEM_FREE(new_vecs_d[i][j], nv, REAL);
825 }
826 for(i = 0; i < 3; i++)
827 if(new_vel[i]) {
828 if(velocity->fe_space->bas_fcts->degree == 0)
829 MEM_FREE(new_vel[i], ne, REAL);
830 else
831 MEM_FREE(new_vel[i], nv, REAL);
832 }
833
834 free_dof_int_vec(dof_vert_ind);
835
836 if(kill_fe_space)
837 free_fe_space(max_fe_space);
838
839 return true;
840 }
841
842
843 /***************************************************************************/
844 /* User interface */
845 /***************************************************************************/
846
847 /***************************************************************************/
848 /* write_mesh_gmv(): */
849 /* Write mesh and vectors in GMV ASCII or binary format. */
850 /* We may output up to 250 DOF_REAL_VECs and 250 DOF_REAL_D_VECs into the */
851 /* GMV file. These are passed as lists. The "velocity" argument is treated */
852 /* in a special way; GMV will automatically generate a field storing the */
853 /* velocity magnitudes. Only the mesh and file names are mandatory. See */
854 /* the description of "add_refined_data()" above for an explanation of */
855 /* "use_refined_grid". The entry "sim_time" represents the numerical time */
856 /* for instationary problems. GMV can display this value together with */
857 /* the FEM data. */
858 /* */
859 /* Return value is 0 for OK and 1 for ERROR. */
860 /***************************************************************************/
861
write_mesh_gmv(MESH * mesh,const char * file_name,bool write_ascii,bool use_refined_grid,const int n_drv,DOF_REAL_VEC ** drv_ptr,const int n_drv_d,DOF_REAL_VEC_D ** drv_d_ptr,DOF_REAL_VEC_D * velocity,REAL sim_time)862 bool write_mesh_gmv(MESH *mesh, const char *file_name,
863 bool write_ascii,
864 bool use_refined_grid,
865 const int n_drv,
866 DOF_REAL_VEC **drv_ptr,
867 const int n_drv_d,
868 DOF_REAL_VEC_D **drv_d_ptr,
869 DOF_REAL_VEC_D *velocity,
870 REAL sim_time)
871 {
872 FUNCNAME("write_mesh_gmv");
873 FILE *file = NULL;
874
875 if(!mesh) {
876 ERROR("no mesh - no file created!\n");
877 return(1);
878 }
879
880 if(n_drv < 0 || n_drv > 250) {
881 ERROR("n_drv must be an int between 0 and 250!\n");
882 return(1);
883 }
884
885 if(n_drv_d < 0 || n_drv_d > 250) {
886 ERROR("n_drv_d must be an int between 0 and 250!\n");
887 return(1);
888 }
889
890 if(write_ascii)
891 file = gmv_open_ascii(file_name, mesh, sim_time);
892 else
893 file = gmv_open_bin(file_name, sizeof(int), sizeof(REAL), sim_time);
894
895 if(!file) {
896 ERROR("cannot open file %s\n",file_name);
897 return(1);
898 }
899
900 dof_compress(mesh);
901
902 add_refined_data(file, mesh, use_refined_grid,
903 write_ascii, true, n_drv, drv_ptr,
904 n_drv_d, drv_d_ptr, velocity);
905
906 if(write_ascii)
907 fprintf(file, "endgmv");
908 else {
909 AI_fwrite("endgmv ",sizeof(char),8,file);
910 }
911
912 fclose(file);
913
914 return(0);
915 }
916
917
918 /***************************************************************************/
919 /* write_dof_vec_gmv(): */
920 /* This function is similar to the one above. However, we do not output */
921 /* the mesh node and element data - instead we direct GMV to access this */
922 /* data from the "mesh_file". This file must have been created by a */
923 /* previous call to write_mesh_gmv(). This way we avoid wasting time and */
924 /* disk space when writing full GMV files including redundant mesh */
925 /* information during each step of an instationary simulation, see GMV */
926 /* documentation. */
927 /* */
928 /* PLEASE NOTE: When using the "use_refined_grid" feature, the prior */
929 /* write_mesh_gmv() call must also be given DOF_REAL_[D_]VECS to */
930 /* determine the needed virtual refinement of the output mesh. */
931 /* */
932 /* Return value is 0 for OK and 1 for ERROR. */
933 /***************************************************************************/
934
write_dof_vec_gmv(MESH * mesh,const char * mesh_file,const char * file_name,bool write_ascii,bool use_refined_grid,const int n_drv,DOF_REAL_VEC ** drv_ptr,const int n_drv_d,DOF_REAL_VEC_D ** drv_d_ptr,DOF_REAL_VEC_D * velocity,REAL sim_time)935 bool write_dof_vec_gmv(MESH *mesh,
936 const char *mesh_file,
937 const char *file_name,
938 bool write_ascii,
939 bool use_refined_grid,
940 const int n_drv,
941 DOF_REAL_VEC **drv_ptr,
942 const int n_drv_d,
943 DOF_REAL_VEC_D **drv_d_ptr,
944 DOF_REAL_VEC_D *velocity,
945 REAL sim_time)
946 {
947 FUNCNAME("write_mesh_gmv");
948 FILE *file = NULL;
949
950 if(n_drv < 0 || n_drv > 250) {
951 ERROR("n_drv must be an int between 0 and 250!\n");
952 return(1);
953 }
954
955 if(n_drv_d < 0 || n_drv_d > 250) {
956 ERROR("n_drv_d must be an int between 0 and 250!\n");
957 return(1);
958 }
959
960
961 if(write_ascii)
962 file = gmv_open_ascii(file_name, mesh, sim_time);
963 else
964 file = gmv_open_bin(file_name, sizeof(int), sizeof(REAL), sim_time);
965
966 if(!file) {
967 ERROR("cannot open file %s\n",file_name);
968 return(1);
969 }
970
971 dof_compress(mesh);
972
973 /* At this point we use the "fromfile" keyword to tell GMV about the other */
974 /* file storing the mesh nodes and elements. */
975
976 if(write_ascii) {
977 fprintf(file, "nodev fromfile \"%s\"\n", mesh_file);
978 fprintf(file, "cells fromfile \"%s\"\n", mesh_file);
979
980 /* Material is only written for parametric meshes. */
981 if(mesh->parametric)
982 fprintf(file, "material fromfile \"%s\"\n", mesh_file);
983 }
984 else {
985 char filename[1024];
986
987 TEST_EXIT(strlen(mesh_file) < 1024, "Sorry, the filename is too long, please use less than 1024 characters.\n");
988
989 snprintf(filename, 1024, "\"%s\"", mesh_file);
990
991 AI_fwrite("nodev fromfile",sizeof(char),16,file);
992 AI_fwrite(filename,sizeof(char),strlen(filename),file);
993 AI_fwrite("cells fromfile",sizeof(char),16,file);
994 AI_fwrite(filename,sizeof(char),strlen(filename),file);
995 /* Material is only written for parametric meshes. */
996 if(mesh->parametric) {
997 AI_fwrite("materialfromfile",sizeof(char),16,file);
998 AI_fwrite(filename,sizeof(char),strlen(filename),file);
999 }
1000 }
1001
1002
1003 add_refined_data(file, mesh, use_refined_grid,
1004 write_ascii, false, n_drv, drv_ptr,
1005 n_drv_d, drv_d_ptr, velocity);
1006
1007 if(write_ascii)
1008 fprintf(file, "endgmv");
1009 else {
1010 AI_fwrite("endgmv ",sizeof(char),8,file);
1011 }
1012
1013 fclose(file);
1014
1015 return(0);
1016 }
1017