/[escript]/branches/scons-dev/finley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /branches/scons-dev/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1655 - (show annotations)
Wed Jul 16 05:10:10 2008 UTC (12 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 40500 byte(s)
Compiles and runs under MPI, and without compiler warnings.
Improved compiler options.
saveVTK() issues a warning (not an error) if VTK is unavailable.

1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 /**************************************************************/
17
18 /* writes data and mesh in a vtk file */
19 /* nodal data needs to be given on FINLEY_NODES or FINLEY_REDUCED_NODES */
20
21 /**************************************************************/
22
23
24 #include "Mesh.h"
25 #include "Assemble.h"
26 #include "vtkCellType.h" /* copied from vtk source directory !!! */
27
28 #define LEN_PRINTED_INT_FORMAT (9+1)
29 #define INT_FORMAT "%d "
30 #define INT_NEWLINE_FORMAT "%d\n"
31 #define FLOAT_SCALAR_FORMAT "%12.6e\n"
32 #define FLOAT_VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
33 #define FLOAT_TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
34 #define LEN_PRINTED_FLOAT_SCALAR_FORMAT (12+1)
35 #define LEN_PRINTED_FLOAT_VECTOR_FORMAT (3*(12+1)+1)
36 #define LEN_PRINTED_FLOAT_TENSOR_FORMAT (9*(12+1)+1)
37 #define NEWLINE "\n"
38 #define LEN_TMP_BUFFER LEN_PRINTED_FLOAT_TENSOR_FORMAT+(MAX_numNodes*LEN_PRINTED_INT_FORMAT+1)+2
39 #define NCOMP_MAX 9
40 #define __STRCAT(dest,chunk,dest_in_use) \
41 { \
42 strcpy(&dest[dest_in_use], chunk); \
43 dest_in_use+=strlen(chunk); \
44 }
45
46 void Finley_Mesh_saveVTK(const char * filename_p,
47 Finley_Mesh *mesh_p,
48 const dim_t num_data,
49 char* *names_p,
50 escriptDataC* *data_pp)
51 {
52 #ifdef USE_VTK
53 char error_msg[LenErrorMsg_MAX], *txt_buffer=NULL, tmp_buffer[LEN_TMP_BUFFER];
54 double sampleAvg[NCOMP_MAX], *values, rtmp;
55 size_t txt_buffer_in_use;
56 dim_t len_txt_buffer, max_len_names;
57 FILE * fileHandle_p = NULL;
58 int mpi_size, i, j, cellType;
59 dim_t i_data;
60 dim_t nDim, globalNumPoints, numCells, globalNumCells, numVTKNodesPerElement;
61 dim_t myNumPoints, numPointsPerSample, rank, nComp, nCompReqd;
62 dim_t shape, NN, numCellFactor, myNumCells, max_name_len;
63 bool_t *isCellCentered=NULL,write_celldata=FALSE,write_pointdata=FALSE;
64 bool_t set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
65 index_t myFirstNode, myLastNode, *globalNodeIndex, k, *node_index, myFirstCell;
66 #ifdef PASO_MPI
67 int ierr;
68 int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY | MPI_MODE_SEQUENTIAL;
69 MPI_File mpi_fileHandle_p;
70 MPI_Status mpi_status;
71 MPI_Request mpi_req;
72 MPI_Info mpi_info=MPI_INFO_NULL;
73 #endif
74 Paso_MPI_rank my_mpi_rank;
75 int nodetype=FINLEY_NODES;
76 int elementtype=FINLEY_UNKNOWN;
77 char elemTypeStr[32];
78 Finley_NodeMapping *nodeMapping=NULL;
79 Finley_ElementFile* elements=NULL;
80 ElementTypeId TypeId;
81
82
83 /****************************************/
84 /* */
85 /* tags in the vtk file */
86
87 char* tags_header="<?xml version=\"1.0\"?>\n" \
88 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n" \
89 "<UnstructuredGrid>\n" \
90 "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
91 "<Points>\n" \
92 "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
93 char *tag_End_DataArray = "</DataArray>\n";
94 char* tag_End_PointData = "</PointData>\n";
95 char* tag_End_CellData = "</CellData>\n";
96 char *footer = "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n";
97 char* tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
98 char* tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
99 char* tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
100 char* tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
101 char *tags_End_Type_And_Cells = "</DataArray>\n</Cells>\n";
102
103 int VTK_QUADRATIC_HEXAHEDRON_INDEX[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15 };
104 /* if there is no mesh we just return */
105 if (mesh_p==NULL) return;
106
107 my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
108 mpi_size = mesh_p->Nodes->MPIInfo->size;
109 nDim = mesh_p->Nodes->numDim;
110
111 if (! ( (nDim ==2) || (nDim == 3) ) ) {
112 Finley_setError(IO_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
113 return;
114 }
115 /*************************************************************************************/
116
117 /* open the file and check handle */
118
119 if (mpi_size > 1) {
120 #ifdef PASO_MPI
121 /* Collective Call */
122 #ifdef MPIO_HINTS
123 MPI_Info_create(&mpi_info);
124 /* MPI_Info_set(mpi_info, "striping_unit", "424288"); */
125 /* MPI_Info_set(mpi_info, "striping_factor", "16"); */
126 /* MPI_Info_set(mpi_info, "collective_buffering", "true"); */
127 /* MPI_Info_set(mpi_info, "cb_block_size", "131072"); */
128 /* MPI_Info_set(mpi_info, "cb_buffer_size", "1048567"); */
129 /* MPI_Info_set(mpi_info, "cb_nodes", "8"); */
130 /* MPI_Info_set(mpi_info, "access_style", "write_once, sequential"); */
131
132 /*XFS only */
133 /* MPI_Info_set(mpi_info, "direct_write", "true"); */
134 #endif
135 ierr=MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p, amode,mpi_info, &mpi_fileHandle_p);
136 if (ierr != MPI_SUCCESS) {
137 perror(filename_p);
138 sprintf(error_msg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
139 Finley_setError(IO_ERROR,error_msg);
140 } else {
141 MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,MPI_CHAR, MPI_CHAR, "native" , mpi_info);
142 }
143 #endif
144 } else {
145 fileHandle_p = fopen(filename_p, "w");
146 if (fileHandle_p==NULL) {
147 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
148 Finley_setError(IO_ERROR,error_msg);
149 }
150 }
151 if (! Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo) ) return;
152 /*************************************************************************************/
153
154 /* find the mesh type to be written */
155
156 isCellCentered=TMPMEMALLOC(num_data,bool_t);
157 max_len_names=0;
158 if (!Finley_checkPtr(isCellCentered)) {
159 nodetype=FINLEY_UNKNOWN;
160 elementtype=FINLEY_UNKNOWN;
161 for (i_data=0;i_data<num_data;++i_data) {
162 if (! isEmpty(data_pp[i_data])) {
163 switch(getFunctionSpaceType(data_pp[i_data]) ) {
164 case FINLEY_NODES:
165 nodetype = (nodetype == FINLEY_REDUCED_NODES) ? FINLEY_REDUCED_NODES : FINLEY_NODES;
166 isCellCentered[i_data]=FALSE;
167 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
168 elementtype=FINLEY_ELEMENTS;
169 } else {
170 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
171 }
172 break;
173 case FINLEY_REDUCED_NODES:
174 nodetype = FINLEY_REDUCED_NODES;
175 isCellCentered[i_data]=FALSE;
176 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
177 elementtype=FINLEY_ELEMENTS;
178 } else {
179 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
180 }
181 break;
182 case FINLEY_ELEMENTS:
183 case FINLEY_REDUCED_ELEMENTS:
184 isCellCentered[i_data]=TRUE;
185 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_ELEMENTS) {
186 elementtype=FINLEY_ELEMENTS;
187 } else {
188 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
189 }
190 break;
191 case FINLEY_FACE_ELEMENTS:
192 case FINLEY_REDUCED_FACE_ELEMENTS:
193 isCellCentered[i_data]=TRUE;
194 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_FACE_ELEMENTS) {
195 elementtype=FINLEY_FACE_ELEMENTS;
196 } else {
197 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
198 }
199 break;
200 case FINLEY_POINTS:
201 isCellCentered[i_data]=TRUE;
202 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_POINTS) {
203 elementtype=FINLEY_POINTS;
204 } else {
205 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
206 }
207 break;
208 case FINLEY_CONTACT_ELEMENTS_1:
209 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
210 isCellCentered[i_data]=TRUE;
211 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
212 elementtype=FINLEY_CONTACT_ELEMENTS_1;
213 } else {
214 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
215 }
216 break;
217 case FINLEY_CONTACT_ELEMENTS_2:
218 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
219 isCellCentered[i_data]=TRUE;
220 if (elementtype==FINLEY_UNKNOWN || elementtype==FINLEY_CONTACT_ELEMENTS_1) {
221 elementtype=FINLEY_CONTACT_ELEMENTS_1;
222 } else {
223 Finley_setError(TYPE_ERROR,"saveVTK: cannot write given data in single file.");
224 }
225 break;
226 default:
227 sprintf(error_msg,"saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[i_data]));
228 Finley_setError(TYPE_ERROR,error_msg);
229 }
230 if (isCellCentered[i_data]) {
231 write_celldata=TRUE;
232 } else {
233 write_pointdata=TRUE;
234 }
235 max_len_names =MAX(max_len_names,(dim_t)strlen(names_p[i_data]));
236 }
237 }
238 nodetype = (nodetype == FINLEY_UNKNOWN) ? FINLEY_NODES : nodetype;
239 }
240 if (Finley_noError()) {
241
242 /***************************************/
243
244 /* select number of points and the mesh component */
245
246 if (nodetype == FINLEY_REDUCED_NODES) {
247 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
248 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
249 globalNumPoints = Finley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
250 globalNodeIndex= Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
251 } else {
252 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
253 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
254 globalNumPoints = Finley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
255 globalNodeIndex= Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
256 }
257 myNumPoints = myLastNode - myFirstNode;
258 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
259 switch(elementtype) {
260 case FINLEY_ELEMENTS:
261 elements=mesh_p->Elements;
262 break;
263 case FINLEY_FACE_ELEMENTS:
264 elements=mesh_p->FaceElements;
265 break;
266 case FINLEY_POINTS:
267 elements=mesh_p->Points;
268 break;
269 case FINLEY_CONTACT_ELEMENTS_1:
270 elements=mesh_p->ContactElements;
271 break;
272 }
273 if (elements==NULL) {
274 Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
275 } else {
276 /* map finley element type to VTK element type */
277 numCells = elements->numElements;
278 globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
279 myNumCells= Finley_ElementFile_getMyNumElements(elements);
280 myFirstCell= Finley_ElementFile_getFirstElement(elements);
281 NN = elements->numNodes;
282 if (nodetype==FINLEY_REDUCED_NODES) {
283 TypeId = elements->LinearReferenceElement->Type->TypeId;
284 } else {
285 TypeId = elements->ReferenceElement->Type->TypeId;
286 }
287 switch(TypeId) {
288 case Point1:
289 case Line2Face:
290 case Line3Face:
291 case Point1_Contact:
292 case Line2Face_Contact:
293 case Line3Face_Contact:
294 numCellFactor=1;
295 cellType = VTK_VERTEX;
296 numVTKNodesPerElement = 1;
297 strcpy(elemTypeStr, "VTK_VERTEX");
298 break;
299
300 case Line2:
301 case Tri3Face:
302 case Rec4Face:
303 case Line2_Contact:
304 case Tri3_Contact:
305 case Tri3Face_Contact:
306 case Rec4Face_Contact:
307 numCellFactor=1;
308 cellType = VTK_LINE;
309 numVTKNodesPerElement = 2;
310 strcpy(elemTypeStr, "VTK_LINE");
311 break;
312
313 case Tri3:
314 case Tet4Face:
315 case Tet4Face_Contact:
316 numCellFactor=1;
317 cellType = VTK_TRIANGLE;
318 numVTKNodesPerElement = 3;
319 strcpy(elemTypeStr, "VTK_TRIANGLE");
320 break;
321
322 case Rec4:
323 case Hex8Face:
324 case Rec4_Contact:
325 case Hex8Face_Contact:
326 numCellFactor=1;
327 cellType = VTK_QUAD;
328 numVTKNodesPerElement = 4;
329 strcpy(elemTypeStr, "VTK_QUAD");
330 break;
331
332 case Tet4:
333 numCellFactor=1;
334 cellType = VTK_TETRA;
335 numVTKNodesPerElement = 4;
336 strcpy(elemTypeStr, "VTK_TETRA");
337 break;
338
339 case Hex8:
340 numCellFactor=1;
341 cellType = VTK_HEXAHEDRON;
342 numVTKNodesPerElement = 8;
343 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
344 break;
345
346 case Line3:
347 case Tri6Face:
348 case Rec8Face:
349 case Line3_Contact:
350 case Tri6Face_Contact:
351 case Rec8Face_Contact:
352 numCellFactor=1;
353 cellType = VTK_QUADRATIC_EDGE;
354 numVTKNodesPerElement = 3;
355 strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
356 break;
357
358 case Tri6:
359 case Tet10Face:
360 case Tri6_Contact:
361 case Tet10Face_Contact:
362 numCellFactor=1;
363 cellType = VTK_QUADRATIC_TRIANGLE;
364 numVTKNodesPerElement = 6;
365 strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
366 break;
367
368 case Rec8:
369 case Hex20Face:
370 case Rec8_Contact:
371 case Hex20Face_Contact:
372 numCellFactor=1;
373 cellType = VTK_QUADRATIC_QUAD;
374 numVTKNodesPerElement = 8;
375 strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
376 break;
377
378 case Tet10:
379 numCellFactor=1;
380 cellType = VTK_QUADRATIC_TETRA;
381 numVTKNodesPerElement = 10;
382 strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
383 break;
384
385 case Hex20:
386 numCellFactor=1;
387 cellType = VTK_QUADRATIC_HEXAHEDRON;
388 numVTKNodesPerElement = 20;
389 strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
390 break;
391
392 default:
393 sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
394 Finley_setError(VALUE_ERROR,error_msg);
395 }
396 }
397 }
398 /***************************************/
399
400 /***************************************/
401 /* */
402 /* allocate text buffer */
403 /* */
404 max_name_len=0;
405 for (i_data =0 ;i_data<num_data;++i_data) max_name_len=MAX(max_name_len,(dim_t)strlen(names_p[i_data]));
406 len_txt_buffer= strlen(tags_header) + 3 * LEN_PRINTED_INT_FORMAT + (30+3*max_name_len); /* header */
407 if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints * LEN_TMP_BUFFER);
408 if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*(LEN_PRINTED_INT_FORMAT*numVTKNodesPerElement+1));
409 len_txt_buffer=MAX(len_txt_buffer,200+3*max_len_names);
410 len_txt_buffer=MAX(len_txt_buffer, (dim_t)strlen(tag_Float_DataArray) + LEN_PRINTED_INT_FORMAT + max_len_names);
411 if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, numCellFactor*myNumCells*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
412 if (mpi_size > 1) len_txt_buffer=MAX(len_txt_buffer, myNumPoints*LEN_PRINTED_FLOAT_TENSOR_FORMAT);
413 txt_buffer=TMPMEMALLOC(len_txt_buffer+1,char);
414 Finley_checkPtr(txt_buffer);
415
416 if (Finley_noError()) {
417
418 /* select number of points and the mesh component */
419
420 sprintf(txt_buffer,tags_header,globalNumPoints,numCellFactor*globalNumCells,3);
421
422 if (mpi_size > 1) {
423 if ( my_mpi_rank == 0) {
424 #ifdef PASO_MPI
425 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
426 MPI_Wait(&mpi_req,&mpi_status);
427 #endif
428 }
429 } else {
430 fprintf(fileHandle_p,txt_buffer);
431 }
432
433 /* write the nodes */
434
435 if (mpi_size > 1) {
436
437 txt_buffer[0] = '\0';
438 txt_buffer_in_use=0;
439 if (nDim==2) {
440 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
441 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
442 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
443 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
444 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
445 0.);
446 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
447 }
448 }
449 } else {
450 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
451 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
452 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,
453 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
454 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
455 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
456 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
457 }
458 }
459
460 }
461 #ifdef PASO_MPI
462 MPI_File_write_ordered(mpi_fileHandle_p, txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
463 #endif
464 } else {
465 if (nDim==2) {
466 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
467 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
468 fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
469 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
470 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
471 0.);
472 }
473 }
474 } else {
475 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
476 if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
477 fprintf(fileHandle_p,FLOAT_VECTOR_FORMAT,
478 mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
479 mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
480 mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
481 }
482 }
483
484 }
485 }
486
487 /* close the Points and open connectivity */
488
489 if (mpi_size > 1) {
490 if ( my_mpi_rank == 0) {
491 #ifdef PASO_MPI
492 MPI_File_iwrite_shared(mpi_fileHandle_p, tags_End_Points_and_Start_Conn, strlen(tags_End_Points_and_Start_Conn), MPI_CHAR, &mpi_req);
493 MPI_Wait(&mpi_req,&mpi_status);
494 #endif
495 }
496 } else {
497 fprintf(fileHandle_p,tags_End_Points_and_Start_Conn);
498 }
499
500 /* write the cells */
501 if (nodetype == FINLEY_REDUCED_NODES) {
502 node_index=elements->ReferenceElement->Type->linearNodes;
503 } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
504 node_index=VTK_QUADRATIC_HEXAHEDRON_INDEX;
505 } else if (numVTKNodesPerElement!=NN) {
506 node_index=elements->ReferenceElement->Type->geoNodes;
507 } else {
508 node_index=NULL;
509 }
510
511 if ( mpi_size > 1) {
512 txt_buffer[0] = '\0';
513 txt_buffer_in_use=0;
514 if (node_index == NULL) {
515 for (i = 0; i < numCells; i++) {
516 if (elements->Owner[i] == my_mpi_rank) {
517 for (j = 0; j < numVTKNodesPerElement; j++) {
518 sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
519 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
520 }
521 __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
522 }
523 }
524 } else {
525 for (i = 0; i < numCells; i++) {
526 if (elements->Owner[i] == my_mpi_rank) {
527 for (j = 0; j < numVTKNodesPerElement; j++) {
528 sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
529 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
530 }
531 __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
532 }
533 }
534 }
535 #ifdef PASO_MPI
536 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
537 #endif
538 } else {
539 if (node_index == NULL) {
540 for (i = 0; i < numCells; i++) {
541 for (j = 0; j < numVTKNodesPerElement; j++) {
542 fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
543 }
544 fprintf(fileHandle_p,NEWLINE);
545 }
546 } else {
547 for (i = 0; i < numCells; i++) {
548 for (j = 0; j < numVTKNodesPerElement; j++) {
549 fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
550 }
551 fprintf(fileHandle_p,NEWLINE);
552 }
553 }
554
555 }
556
557 /* finalize the connection and start the offset section */
558 if (mpi_size > 1) {
559 if( my_mpi_rank == 0) {
560 #ifdef PASO_MPI
561 MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Conn_and_Start_Offset,strlen(tags_End_Conn_and_Start_Offset),MPI_CHAR,&mpi_req);
562 MPI_Wait(&mpi_req,&mpi_status);
563 #endif
564 }
565 } else {
566 fprintf(fileHandle_p,tags_End_Conn_and_Start_Offset);
567 }
568
569 /* write the offsets */
570
571 if ( mpi_size > 1) {
572 txt_buffer[0] = '\0';
573 txt_buffer_in_use=0;
574 for (i=numVTKNodesPerElement*(myFirstCell*numCellFactor+1); i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
575 sprintf(tmp_buffer, INT_NEWLINE_FORMAT, i);
576 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
577 }
578 #ifdef PASO_MPI
579 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
580 #endif
581 } else {
582 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
583 fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
584 }
585
586 }
587 /* finalize the offset section and start the type section */
588 if ( mpi_size > 1) {
589 if ( my_mpi_rank == 0) {
590 #ifdef PASO_MPI
591 MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Offset_and_Start_Type,strlen(tags_End_Offset_and_Start_Type),MPI_CHAR,&mpi_req);
592 MPI_Wait(&mpi_req,&mpi_status);
593 #endif
594 }
595 } else {
596 fprintf(fileHandle_p,tags_End_Offset_and_Start_Type);
597 }
598 /* write element type */
599 sprintf(tmp_buffer, INT_NEWLINE_FORMAT, cellType);
600 if ( mpi_size > 1) {
601 txt_buffer[0] = '\0';
602 txt_buffer_in_use=0;
603 for (i=0; i<numCells*numCellFactor; i++) __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
604 #ifdef PASO_MPI
605 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
606 #endif
607 } else {
608 for (i=0; i<numCells*numCellFactor; i++) fprintf(fileHandle_p, tmp_buffer);
609 }
610 /* finalize cell information */
611 if ( mpi_size > 1) {
612 if ( my_mpi_rank == 0) {
613 #ifdef PASO_MPI
614 MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Type_And_Cells,strlen(tags_End_Type_And_Cells),MPI_CHAR,&mpi_req);
615 MPI_Wait(&mpi_req,&mpi_status);
616 #endif
617 }
618 } else {
619 fprintf(fileHandle_p,tags_End_Type_And_Cells);
620 }
621 }
622
623 /* Write cell data */
624 if (write_celldata && Finley_noError()) {
625 /* mark the active data arrays */
626 txt_buffer[0] = '\0';
627 set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
628 strcat(txt_buffer, "<CellData");
629 for (i_data =0 ;i_data<num_data;++i_data) {
630 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
631 /* if the rank == 0: --> scalar data */
632 /* if the rank == 1: --> vector data */
633 /* if the rank == 2: --> tensor data */
634
635 switch(getDataPointRank(data_pp[i_data])) {
636 case 0:
637 if (! set_scalar) {
638 strcat(txt_buffer," Scalars=\"");
639 strcat(txt_buffer,names_p[i_data]);
640 strcat(txt_buffer,"\"");
641 set_scalar=TRUE;
642 }
643 break;
644 case 1:
645 if (! set_vector) {
646 strcat(txt_buffer," Vectors=\"");
647 strcat(txt_buffer,names_p[i_data]);
648 strcat(txt_buffer,"\"");
649 set_vector=TRUE;
650 }
651 break;
652 case 2:
653 if (! set_tensor) {
654 strcat(txt_buffer," Tensors=\"");
655 strcat(txt_buffer,names_p[i_data]);
656 strcat(txt_buffer,"\"");
657 set_tensor=TRUE;
658 }
659 break;
660 default:
661 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
662 Finley_setError(VALUE_ERROR,error_msg);
663 return;
664 }
665 }
666 }
667 strcat(txt_buffer, ">\n");
668 if ( mpi_size > 1) {
669 if ( my_mpi_rank == 0) {
670 #ifdef PASO_MPI
671 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
672 MPI_Wait(&mpi_req,&mpi_status);
673 #endif
674 }
675 } else {
676 fprintf(fileHandle_p,txt_buffer);
677 }
678 /* write the arrays */
679 for (i_data =0 ;i_data<num_data;++i_data) {
680 if (! isEmpty(data_pp[i_data]) && isCellCentered[i_data]) {
681 txt_buffer[0] = '\0';
682 txt_buffer_in_use=0;
683 numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
684 rank = getDataPointRank(data_pp[i_data]);
685 nComp = getDataPointSize(data_pp[i_data]);
686 nCompReqd=1; /* the number of components mpi_required by vtk */
687 shape=0;
688 if (rank == 0) {
689 nCompReqd = 1;
690 } else if (rank == 1) {
691 shape=getDataPointShape(data_pp[i_data], 0);
692 if (shape>3) {
693 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
694 }
695 nCompReqd = 3;
696 } else {
697 shape=getDataPointShape(data_pp[i_data], 0);
698 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
699 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
700 }
701 nCompReqd = 9;
702 }
703 if (Finley_noError()) {
704 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
705 if ( mpi_size > 1) {
706 if ( my_mpi_rank == 0) {
707 #ifdef PASO_MPI
708 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
709 MPI_Wait(&mpi_req,&mpi_status);
710 #endif
711 }
712 } else {
713 fprintf(fileHandle_p,txt_buffer);
714 }
715 for (i=0; i<numCells; i++) {
716 if (elements->Owner[i] == my_mpi_rank) {
717 values = getSampleData(data_pp[i_data], i);
718 /* averaging over the number of points in the sample */
719 for (k=0; k<MIN(nComp,NCOMP_MAX); k++) {
720 if (isExpanded(data_pp[i_data])) {
721 rtmp = 0.;
722 for (j=0; j<numPointsPerSample; j++) rtmp += values[INDEX2(k,j,nComp)];
723 sampleAvg[k] = rtmp/numPointsPerSample;
724 } else {
725 sampleAvg[k] = values[k];
726 }
727 }
728 /* if the number of mpi_required components is more than the number
729 * of actual components, pad with zeros
730 */
731 /* probably only need to get shape of first element */
732 /* write the data different ways for scalar, vector and tensor */
733 if (nCompReqd == 1) {
734 sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,sampleAvg[0]);
735 } else if (nCompReqd == 3) {
736 if (shape==1) {
737 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],0.,0.);
738 } else if (shape==2) {
739 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],0.);
740 } else if (shape==3) {
741 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2]);
742 }
743 } else if (nCompReqd == 9) {
744 if (shape==1) {
745 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],0.,0.,
746 0.,0.,0.,
747 0.,0.,0.);
748 } else if (shape==2) {
749 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],sampleAvg[1],0.,
750 sampleAvg[2],sampleAvg[3],0.,
751 0.,0.,0.);
752 } else if (shape==3) {
753 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],sampleAvg[1],sampleAvg[2],
754 sampleAvg[3],sampleAvg[4],sampleAvg[5],
755 sampleAvg[6],sampleAvg[7],sampleAvg[8]);
756 }
757 }
758 if ( mpi_size > 1) {
759 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
760 } else {
761 fprintf(fileHandle_p,tmp_buffer);
762 }
763 }
764 }
765 if ( mpi_size > 1) {
766 #ifdef PASO_MPI
767 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
768 #endif
769 if ( my_mpi_rank == 0) {
770 #ifdef PASO_MPI
771 MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
772 MPI_Wait(&mpi_req,&mpi_status);
773 #endif
774 }
775 } else {
776 fprintf(fileHandle_p,tag_End_DataArray);
777 }
778 }
779 }
780 }
781 if ( mpi_size > 1) {
782 if ( my_mpi_rank == 0) {
783 #ifdef PASO_MPI
784 MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_CellData),MPI_CHAR,&mpi_req);
785 MPI_Wait(&mpi_req,&mpi_status);
786 #endif
787 }
788 } else {
789 fprintf(fileHandle_p,tag_End_CellData);
790 }
791 }
792 /* point data */
793 if (write_pointdata && Finley_noError()) {
794 /* mark the active data arrays */
795 set_scalar=FALSE,set_vector=FALSE, set_tensor=FALSE;
796 txt_buffer[0] = '\0';
797 strcat(txt_buffer, "<PointData");
798 for (i_data =0 ;i_data<num_data;++i_data) {
799 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
800 /* if the rank == 0: --> scalar data */
801 /* if the rank == 1: --> vector data */
802 /* if the rank == 2: --> tensor data */
803
804 switch(getDataPointRank(data_pp[i_data])) {
805 case 0:
806 if (! set_scalar) {
807 strcat(txt_buffer," Scalars=\"");
808 strcat(txt_buffer,names_p[i_data]);
809 strcat(txt_buffer,"\"");
810 set_scalar=TRUE;
811 }
812 break;
813 case 1:
814 if (! set_vector) {
815 strcat(txt_buffer," Vectors=\"");
816 strcat(txt_buffer,names_p[i_data]);
817 strcat(txt_buffer,"\"");
818 set_vector=TRUE;
819 }
820 break;
821 case 2:
822 if (! set_tensor) {
823 strcat(txt_buffer," Tensors=\"");
824 strcat(txt_buffer,names_p[i_data]);
825 strcat(txt_buffer,"\"");
826 set_tensor=TRUE;
827 }
828 break;
829 default:
830 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
831 Finley_setError(VALUE_ERROR,error_msg);
832 return;
833 }
834 }
835 }
836 strcat(txt_buffer, ">\n");
837 if ( mpi_size > 1) {
838 if ( my_mpi_rank == 0) {
839 #ifdef PASO_MPI
840 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
841 MPI_Wait(&mpi_req,&mpi_status);
842 #endif
843 }
844 } else {
845 fprintf(fileHandle_p,txt_buffer);
846 }
847 /* write the arrays */
848 for (i_data =0 ;i_data<num_data;++i_data) {
849 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
850 txt_buffer[0] = '\0';
851 txt_buffer_in_use=0;
852 numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
853 rank = getDataPointRank(data_pp[i_data]);
854 nComp = getDataPointSize(data_pp[i_data]);
855 if (getFunctionSpaceType(data_pp[i_data]) == FINLEY_REDUCED_NODES) {
856 nodeMapping=mesh_p->Nodes->reducedNodesMapping;
857 } else {
858 nodeMapping=mesh_p->Nodes->nodesMapping;
859 }
860 nCompReqd=1; /* the number of components mpi_required by vtk */
861 shape=0;
862 if (rank == 0) {
863 nCompReqd = 1;
864 } else if (rank == 1) {
865 shape=getDataPointShape(data_pp[i_data], 0);
866 if (shape>3) {
867 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
868 }
869 nCompReqd = 3;
870 } else {
871 shape=getDataPointShape(data_pp[i_data], 0);
872 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
873 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
874 }
875 nCompReqd = 9;
876 }
877 if (Finley_noError()) {
878 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
879 if ( mpi_size > 1) {
880 if ( my_mpi_rank == 0) {
881 #ifdef PASO_MPI
882 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
883 MPI_Wait(&mpi_req,&mpi_status);
884 #endif
885 }
886 } else {
887 fprintf(fileHandle_p,txt_buffer);
888 }
889 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
890 k=globalNodeIndex[i];
891 if ( (myFirstNode <= k) && (k < myLastNode) ) {
892 values = getSampleData(data_pp[i_data], nodeMapping->target[i]);
893 /* if the number of mpi_required components is more than the number
894 * of actual components, pad with zeros
895 */
896 /* probably only need to get shape of first element */
897 /* write the data different ways for scalar, vector and tensor */
898 if (nCompReqd == 1) {
899 sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,values[0]);
900 } else if (nCompReqd == 3) {
901 if (shape==1) {
902 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],0.,0.);
903 } else if (shape==2) {
904 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],0.);
905 } else if (shape==3) {
906 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],values[2]);
907 }
908 } else if (nCompReqd == 9) {
909 if (shape==1) {
910 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],0.,0.,
911 0.,0.,0.,
912 0.,0.,0.);
913 } else if (shape==2) {
914 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],0.,
915 values[2],values[3],0.,
916 0.,0.,0.);
917 } else if (shape==3) {
918 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],values[2],
919 values[3],values[4],values[5],
920 values[6],values[7],values[8]);
921 }
922 }
923 if ( mpi_size > 1) {
924 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
925 } else {
926 fprintf(fileHandle_p,tmp_buffer);
927 }
928 }
929 }
930 if ( mpi_size > 1) {
931 #ifdef PASO_MPI
932 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
933 #endif
934 if ( my_mpi_rank == 0) {
935 #ifdef PASO_MPI
936 MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
937 MPI_Wait(&mpi_req,&mpi_status);
938 #endif
939 }
940 } else {
941 fprintf(fileHandle_p,tag_End_DataArray);
942 }
943 }
944 }
945 }
946 if ( mpi_size > 1) {
947 if ( my_mpi_rank == 0) {
948 #ifdef PASO_MPI
949 MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_PointData),MPI_CHAR,&mpi_req);
950 MPI_Wait(&mpi_req,&mpi_status);
951 #endif
952 }
953 } else {
954 fprintf(fileHandle_p,tag_End_PointData);
955 }
956 }
957 if (Finley_noError()) {
958 if ( mpi_size > 1) {
959 if ( my_mpi_rank == 0) {
960 #ifdef PASO_MPI
961 MPI_File_iwrite_shared(mpi_fileHandle_p,footer,strlen(footer),MPI_CHAR,&mpi_req);
962 MPI_Wait(&mpi_req,&mpi_status);
963 #ifdef MPIO_HINTS
964 MPI_Info_free(&mpi_info);
965 #undef MPIO_HINTS
966 #endif
967 #endif
968 }
969 #ifdef PASO_MPI
970 MPI_File_close(&mpi_fileHandle_p);
971 #endif
972 } else {
973 fprintf(fileHandle_p,footer);
974 fclose(fileHandle_p);
975 }
976 }
977 TMPMEMFREE(isCellCentered);
978 TMPMEMFREE(txt_buffer);
979 return;
980 #else
981 /* Don't kill the job if saveVTK() doesn't work */
982 fprintf(stderr, "\n\nsaveVTK: VTK is not available\n\n\n");
983 #endif
984 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26