/[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 1562 - (show annotations)
Wed May 21 13:04:40 2008 UTC (11 years, 5 months ago) by gross
File MIME type: text/plain
File size: 40470 byte(s)
The algebraic upwinding with MPI. The case of boundary constraint needs still some attention. 


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26