/[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 1677 - (show annotations)
Mon Jul 28 06:36:19 2008 UTC (11 years ago) by gross
File MIME type: text/plain
File size: 40603 byte(s)
bug in the MPI version of saveVTK fixed
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<=(myFirstCell+myNumCells)*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