/[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 1274 - (show annotations)
Fri Aug 24 03:09:22 2007 UTC (12 years, 9 months ago) by gross
File MIME type: text/plain
File size: 40931 byte(s)
the vtk tests pass now. dx is still failing
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 "%12.6e\n"
37 #define FLOAT_VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
38 #define FLOAT_TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
39 #define LEN_PRINTED_FLOAT_SCALAR_FORMAT (12+1)
40 #define LEN_PRINTED_FLOAT_VECTOR_FORMAT (3*(12+1)+1)
41 #define LEN_PRINTED_FLOAT_TENSOR_FORMAT (9*(12+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, max_name_len;
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 } else {
253 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
254 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
255 globalNumPoints = Finley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
256 globalNodeIndex= Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
257 }
258 myNumPoints = myLastNode - myFirstNode;
259 if (elementtype==FINLEY_UNKNOWN) elementtype=FINLEY_ELEMENTS;
260 switch(elementtype) {
261 case FINLEY_ELEMENTS:
262 elements=mesh_p->Elements;
263 break;
264 case FINLEY_FACE_ELEMENTS:
265 elements=mesh_p->FaceElements;
266 break;
267 case FINLEY_POINTS:
268 elements=mesh_p->Points;
269 break;
270 case FINLEY_CONTACT_ELEMENTS_1:
271 elements=mesh_p->ContactElements;
272 break;
273 }
274 if (elements==NULL) {
275 Finley_setError(SYSTEM_ERROR,"saveVTK: undefined element file");
276 } else {
277 /* map finley element type to VTK element type */
278 numCells = elements->numElements;
279 globalNumCells = Finley_ElementFile_getGlobalNumElements(elements);
280 myNumCells= Finley_ElementFile_getMyNumElements(elements);
281 myFirstCell= Finley_ElementFile_getFirstElement(elements);
282 NN = elements->numNodes;
283 if (nodetype==FINLEY_REDUCED_NODES || nodetype==FINLEY_REDUCED_NODES) {
284 TypeId = elements->LinearReferenceElement->Type->TypeId;
285 } else {
286 TypeId = elements->ReferenceElement->Type->TypeId;
287 }
288 switch(TypeId) {
289 case Point1:
290 case Line2Face:
291 case Line3Face:
292 case Point1_Contact:
293 case Line2Face_Contact:
294 case Line3Face_Contact:
295 numCellFactor=1;
296 cellType = VTK_VERTEX;
297 numVTKNodesPerElement = 1;
298 strcpy(elemTypeStr, "VTK_VERTEX");
299 break;
300
301 case Line2:
302 case Tri3Face:
303 case Rec4Face:
304 case Line2_Contact:
305 case Tri3_Contact:
306 case Tri3Face_Contact:
307 case Rec4Face_Contact:
308 numCellFactor=1;
309 cellType = VTK_LINE;
310 numVTKNodesPerElement = 2;
311 strcpy(elemTypeStr, "VTK_LINE");
312 break;
313
314 case Tri3:
315 case Tet4Face:
316 case Tet4Face_Contact:
317 numCellFactor=1;
318 cellType = VTK_TRIANGLE;
319 numVTKNodesPerElement = 3;
320 strcpy(elemTypeStr, "VTK_TRIANGLE");
321 break;
322
323 case Rec4:
324 case Hex8Face:
325 case Rec4_Contact:
326 case Hex8Face_Contact:
327 numCellFactor=1;
328 cellType = VTK_QUAD;
329 numVTKNodesPerElement = 4;
330 strcpy(elemTypeStr, "VTK_QUAD");
331 break;
332
333 case Tet4:
334 numCellFactor=1;
335 cellType = VTK_TETRA;
336 numVTKNodesPerElement = 4;
337 strcpy(elemTypeStr, "VTK_TETRA");
338 break;
339
340 case Hex8:
341 numCellFactor=1;
342 cellType = VTK_HEXAHEDRON;
343 numVTKNodesPerElement = 8;
344 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
345 break;
346
347 case Line3:
348 case Tri6Face:
349 case Rec8Face:
350 case Line3_Contact:
351 case Tri6Face_Contact:
352 case Rec8Face_Contact:
353 numCellFactor=1;
354 cellType = VTK_QUADRATIC_EDGE;
355 numVTKNodesPerElement = 3;
356 strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
357 break;
358
359 case Tri6:
360 case Tet10Face:
361 case Tri6_Contact:
362 case Tet10Face_Contact:
363 numCellFactor=1;
364 cellType = VTK_QUADRATIC_TRIANGLE;
365 numVTKNodesPerElement = 6;
366 strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
367 break;
368
369 case Rec8:
370 case Hex20Face:
371 case Rec8_Contact:
372 case Hex20Face_Contact:
373 numCellFactor=1;
374 cellType = VTK_QUADRATIC_QUAD;
375 numVTKNodesPerElement = 8;
376 strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
377 break;
378
379 case Tet10:
380 numCellFactor=1;
381 cellType = VTK_QUADRATIC_TETRA;
382 numVTKNodesPerElement = 10;
383 strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
384 break;
385
386 case Hex20:
387 numCellFactor=1;
388 cellType = VTK_QUADRATIC_HEXAHEDRON;
389 numVTKNodesPerElement = 20;
390 strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
391 break;
392
393 default:
394 sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
395 Finley_setError(VALUE_ERROR,error_msg);
396 }
397 }
398 }
399 /***************************************/
400
401 /***************************************/
402 /* */
403 /* allocate text buffer */
404 /* */
405 max_name_len=0;
406 for (i_data =0 ;i_data<num_data;++i_data) max_name_len=MAX(max_name_len,strlen(names_p[i_data]));
407 len_txt_buffer= strlen(tags_header) + 3 * LEN_PRINTED_INT_FORMAT + (30+3*max_name_len); /* 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 }
497 } else {
498 fprintf(fileHandle_p,tags_End_Points_and_Start_Conn);
499 }
500
501 /* write the cells */
502 if (nodetype == FINLEY_REDUCED_NODES) {
503 node_index=elements->ReferenceElement->Type->linearNodes;
504 } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
505 node_index=VTK_QUADRATIC_HEXAHEDRON_INDEX;
506 } else if (numVTKNodesPerElement!=NN) {
507 node_index=elements->ReferenceElement->Type->geoNodes;
508 } else {
509 node_index=NULL;
510 }
511
512 if ( mpi_size > 1) {
513 txt_buffer[0] = '\0';
514 txt_buffer_in_use=0;
515 if (node_index == NULL) {
516 for (i = 0; i < numCells; i++) {
517 if (elements->Owner[i] == my_mpi_rank) {
518 for (j = 0; j < numVTKNodesPerElement; j++) {
519 sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
520 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
521 }
522 __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
523 }
524 }
525 } else {
526 for (i = 0; i < numCells; i++) {
527 if (elements->Owner[i] == my_mpi_rank) {
528 for (j = 0; j < numVTKNodesPerElement; j++) {
529 sprintf(tmp_buffer,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
530 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use)
531 }
532 __STRCAT(txt_buffer,NEWLINE,txt_buffer_in_use)
533 }
534 }
535 }
536 #ifdef PASO_MPI
537 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
538 #endif
539 } else {
540 if (node_index == NULL) {
541 for (i = 0; i < numCells; i++) {
542 for (j = 0; j < numVTKNodesPerElement; j++) {
543 fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
544 }
545 fprintf(fileHandle_p,NEWLINE);
546 }
547 } else {
548 for (i = 0; i < numCells; i++) {
549 for (j = 0; j < numVTKNodesPerElement; j++) {
550 fprintf(fileHandle_p,INT_FORMAT,globalNodeIndex[elements->Nodes[INDEX2(node_index[j], i, NN)]]);
551 }
552 fprintf(fileHandle_p,NEWLINE);
553 }
554 }
555
556 }
557
558 /* finalize the connection and start the offset section */
559 if (mpi_size > 1) {
560 if( my_mpi_rank == 0) {
561 #ifdef PASO_MPI
562 MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Conn_and_Start_Offset,strlen(tags_End_Conn_and_Start_Offset),MPI_CHAR,&mpi_req);
563 MPI_Wait(&mpi_req,&mpi_status);
564 #endif
565 }
566 } else {
567 fprintf(fileHandle_p,tags_End_Conn_and_Start_Offset);
568 }
569
570 /* write the offsets */
571
572 if ( mpi_size > 1) {
573 txt_buffer[0] = '\0';
574 txt_buffer_in_use=0;
575 for (i=numVTKNodesPerElement*(myFirstCell*numCellFactor+1); i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
576 sprintf(tmp_buffer, INT_NEWLINE_FORMAT, i);
577 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
578 }
579 #ifdef PASO_MPI
580 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
581 #endif
582 } else {
583 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement*numCellFactor; i+=numVTKNodesPerElement) {
584 fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
585 }
586
587 }
588 /* finalize the offset section and start the type section */
589 if ( mpi_size > 1) {
590 if ( my_mpi_rank == 0) {
591 #ifdef PASO_MPI
592 MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Offset_and_Start_Type,strlen(tags_End_Offset_and_Start_Type),MPI_CHAR,&mpi_req);
593 MPI_Wait(&mpi_req,&mpi_status);
594 #endif
595 }
596 } else {
597 fprintf(fileHandle_p,tags_End_Offset_and_Start_Type);
598 }
599
600
601 /* write element type */
602 sprintf(tmp_buffer, INT_NEWLINE_FORMAT, cellType);
603 if ( mpi_size > 1) {
604 txt_buffer[0] = '\0';
605 txt_buffer_in_use=0;
606 for (i=0; i<numCells*numCellFactor; i++) __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
607 #ifdef PASO_MPI
608 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
609 #endif
610 } else {
611 for (i=0; i<numCells*numCellFactor; i++) fprintf(fileHandle_p, tmp_buffer);
612 }
613 /* finalize cell information */
614 if ( mpi_size > 1) {
615 if ( my_mpi_rank == 0) {
616 #ifdef PASO_MPI
617 MPI_File_iwrite_shared(mpi_fileHandle_p,tags_End_Type_And_Cells,strlen(tags_End_Type_And_Cells),MPI_CHAR,&mpi_req);
618 MPI_Wait(&mpi_req,&mpi_status);
619 #endif
620 }
621 } else {
622 fprintf(fileHandle_p,tags_End_Type_And_Cells);
623 }
624 }
625
626 /* Write cell data */
627 if (write_celldata && Finley_noError()) {
628 /* mark the active data arrays */
629 txt_buffer[0] = '\0';
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[1],sampleAvg[2]);
745 }
746 } else if (nCompReqd == 9) {
747 if (shape==1) {
748 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,sampleAvg[0],0.,0.,
749 0.,0.,0.,
750 0.,0.,0.);
751 } else if (shape==2) {
752 sprintf(tmp_buffer,FLOAT_TENSOR_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_TENSOR_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 txt_buffer[0] = '\0';
800 strcat(txt_buffer, "<PointData");
801 for (i_data =0 ;i_data<num_data;++i_data) {
802 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
803 /* if the rank == 0: --> scalar data */
804 /* if the rank == 1: --> vector data */
805 /* if the rank == 2: --> tensor data */
806
807 switch(getDataPointRank(data_pp[i_data])) {
808 case 0:
809 if (! set_scalar) {
810 strcat(txt_buffer," Scalars=\"");
811 strcat(txt_buffer,names_p[i_data]);
812 strcat(txt_buffer,"\"");
813 set_scalar=TRUE;
814 }
815 break;
816 case 1:
817 if (! set_vector) {
818 strcat(txt_buffer," Vectors=\"");
819 strcat(txt_buffer,names_p[i_data]);
820 strcat(txt_buffer,"\"");
821 set_vector=TRUE;
822 }
823 break;
824 case 2:
825 if (! set_tensor) {
826 strcat(txt_buffer," Tensors=\"");
827 strcat(txt_buffer,names_p[i_data]);
828 strcat(txt_buffer,"\"");
829 set_tensor=TRUE;
830 }
831 break;
832 default:
833 sprintf(error_msg, "saveVTK: data %s: Vtk can't handle objects with rank greater than 2.",names_p[i_data]);
834 Finley_setError(VALUE_ERROR,error_msg);
835 return;
836 }
837 }
838 }
839 strcat(txt_buffer, ">\n");
840 if ( mpi_size > 1) {
841 if ( my_mpi_rank == 0) {
842 #ifdef PASO_MPI
843 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
844 MPI_Wait(&mpi_req,&mpi_status);
845 #endif
846 }
847 } else {
848 fprintf(fileHandle_p,txt_buffer);
849 }
850 /* write the arrays */
851 for (i_data =0 ;i_data<num_data;++i_data) {
852 if (! isEmpty(data_pp[i_data]) && !isCellCentered[i_data]) {
853 txt_buffer[0] = '\0';
854 txt_buffer_in_use=0;
855 numPointsPerSample=getNumDataPointsPerSample(data_pp[i_data]);
856 rank = getDataPointRank(data_pp[i_data]);
857 nComp = getDataPointSize(data_pp[i_data]);
858 if (getFunctionSpaceType(data_pp[i_data]) == FINLEY_REDUCED_NODES) {
859 nodeMapping=mesh_p->Nodes->reducedNodesMapping;
860 } else {
861 nodeMapping=mesh_p->Nodes->nodesMapping;
862 }
863 nCompReqd=1; /* the number of components mpi_required by vtk */
864 shape=0;
865 if (rank == 0) {
866 nCompReqd = 1;
867 } else if (rank == 1) {
868 shape=getDataPointShape(data_pp[i_data], 0);
869 if (shape>3) {
870 Finley_setError(VALUE_ERROR, "saveVTK: rank 1 object must have less then 4 components");
871 }
872 nCompReqd = 3;
873 } else {
874 shape=getDataPointShape(data_pp[i_data], 0);
875 if (shape>3 || shape != getDataPointShape(data_pp[i_data], 1)) {
876 Finley_setError(VALUE_ERROR, "saveVTK: rank 2 object must have less then 4x4 components and must have a square shape");
877 }
878 nCompReqd = 9;
879 }
880 if (Finley_noError()) {
881 sprintf(txt_buffer,tag_Float_DataArray,names_p[i_data], nCompReqd);
882 if ( mpi_size > 1) {
883 if ( my_mpi_rank == 0) {
884 #ifdef PASO_MPI
885 MPI_File_iwrite_shared(mpi_fileHandle_p,txt_buffer,strlen(txt_buffer),MPI_CHAR,&mpi_req);
886 MPI_Wait(&mpi_req,&mpi_status);
887 #endif
888 }
889 } else {
890 fprintf(fileHandle_p,txt_buffer);
891 }
892 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
893 k=globalNodeIndex[i];
894 if ( (myFirstNode <= k) && (k < myLastNode) ) {
895 values = getSampleData(data_pp[i_data], nodeMapping->target[i]);
896 /* if the number of mpi_required components is more than the number
897 * of actual components, pad with zeros
898 */
899 /* probably only need to get shape of first element */
900 /* write the data different ways for scalar, vector and tensor */
901 if (nCompReqd == 1) {
902 sprintf(tmp_buffer,FLOAT_SCALAR_FORMAT,values[0]);
903 } else if (nCompReqd == 3) {
904 if (shape==1) {
905 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],0.,0.);
906 } else if (shape==2) {
907 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],0.);
908 } else if (shape==3) {
909 sprintf(tmp_buffer,FLOAT_VECTOR_FORMAT,values[0],values[1],values[2]);
910 }
911 } else if (nCompReqd == 9) {
912 if (shape==1) {
913 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],0.,0.,
914 0.,0.,0.,
915 0.,0.,0.);
916 } else if (shape==2) {
917 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],0.,
918 values[2],values[3],0.,
919 0.,0.,0.);
920 } else if (shape==3) {
921 sprintf(tmp_buffer,FLOAT_TENSOR_FORMAT,values[0],values[1],values[2],
922 values[3],values[4],values[5],
923 values[6],values[7],values[8]);
924 }
925 }
926 if ( mpi_size > 1) {
927 __STRCAT(txt_buffer,tmp_buffer,txt_buffer_in_use);
928 } else {
929 fprintf(fileHandle_p,tmp_buffer);
930 }
931 }
932 }
933 if ( mpi_size > 1) {
934 #ifdef PASO_MPI
935 MPI_File_write_ordered(mpi_fileHandle_p,txt_buffer,txt_buffer_in_use, MPI_CHAR, &mpi_status);
936 #endif
937 if ( my_mpi_rank == 0) {
938 #ifdef PASO_MPI
939 MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_DataArray,strlen(tag_End_DataArray),MPI_CHAR,&mpi_req);
940 MPI_Wait(&mpi_req,&mpi_status);
941 #endif
942 }
943 } else {
944 fprintf(fileHandle_p,tag_End_DataArray);
945 }
946 }
947 }
948 }
949 if ( mpi_size > 1) {
950 if ( my_mpi_rank == 0) {
951 #ifdef PASO_MPI
952 MPI_File_iwrite_shared(mpi_fileHandle_p,tag_End_CellData,strlen(tag_End_PointData),MPI_CHAR,&mpi_req);
953 MPI_Wait(&mpi_req,&mpi_status);
954 #endif
955 }
956 } else {
957 fprintf(fileHandle_p,tag_End_PointData);
958 }
959 }
960 if (Finley_noError()) {
961 if ( mpi_size > 1) {
962 if ( my_mpi_rank == 0) {
963 #ifdef PASO_MPI
964 MPI_File_iwrite_shared(mpi_fileHandle_p,footer,strlen(footer),MPI_CHAR,&mpi_req);
965 MPI_Wait(&mpi_req,&mpi_status);
966 #ifdef MPIO_HINTS
967 MPI_Info_free(&mpi_info);
968 #undef MPIO_HINTS
969 #endif
970 MPI_File_close(&mpi_fileHandle_p);
971 #endif
972 }
973 } else {
974 fprintf(fileHandle_p,footer);
975 fclose(fileHandle_p);
976 }
977 }
978 TMPMEMFREE(isCellCentered);
979 TMPMEMFREE(txt_buffer);
980 return;
981 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26