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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26