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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26