/[escript]/trunk/finley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1669 - (show annotations)
Thu Jul 24 01:10:04 2008 UTC (11 years, 2 months ago) by gross
File MIME type: text/plain
File size: 40585 byte(s)
A problem with VTK writer and MPI is fixed: Apparently MPI_file_open does not "delete" the file  which has the effect that 
if less date are written into the file as the file contained when opened bits of the previous containt remains in the file. This problem is fixed
by deleting the file if it exists before open it with MPI. The additional function Paso_existFile needed to be added as there is no standart C 
function to test the exists of a file.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26