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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26