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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1675 - (show annotations)
Fri Jul 25 05:48:29 2008 UTC (14 years, 8 months ago) by ksteube
File MIME type: text/plain
File size: 40736 byte(s)
Tweaked the default settings for various options to make configuration more automatic.
Fixed configuration for MKL.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26