/[escript]/trunk-mpi-branch/finley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /trunk-mpi-branch/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1223 - (show annotations)
Fri Aug 3 02:40:39 2007 UTC (11 years, 8 months ago) by gross
File MIME type: text/plain
File size: 40589 byte(s)
first attemt towards an improved MPI version.  

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