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