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