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

Annotation of /trunk/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1388 - (hide annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 10 months ago) by trankine
File MIME type: text/plain
File size: 40308 byte(s)
And get the *(&(*&(* name right
1 jgs 110
2 ksteube 1312 /* $Id$ */
3 jgs 150
4 ksteube 1312 /*******************************************************
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 jgs 110 /**************************************************************/
17    
18     /* writes data and mesh in a vtk file */
19 ksteube 1312 /* nodal data needs to be given on FINLEY_NODES or FINLEY_REDUCED_NODES */
20 jgs 110
21     /**************************************************************/
22    
23 dhawcroft 793
24 jgs 110 #include "Mesh.h"
25 gross 1062 #include "Assemble.h"
26 jgs 113 #include "vtkCellType.h" /* copied from vtk source directory !!! */
27 jgs 110
28 ksteube 1312 #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 dhawcroft 793 }
45    
46 ksteube 1312 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 dhawcroft 793 {
52 ksteube 1312 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 dhawcroft 793 int amode = MPI_MODE_CREATE | MPI_MODE_WRONLY | MPI_MODE_SEQUENTIAL;
64 ksteube 1312 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 dhawcroft 793 int elementtype=FINLEY_UNKNOWN;
72     char elemTypeStr[32];
73 ksteube 1312 Finley_NodeMapping *nodeMapping=NULL;
74 dhawcroft 793 Finley_ElementFile* elements=NULL;
75 ksteube 1312 ElementTypeId TypeId;
76    
77    
78     /****************************************/
79     /* */
80     /* tags in the vtk file */
81 dhawcroft 793
82 ksteube 1312 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 dhawcroft 793
98 ksteube 1312 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 dhawcroft 793
102 ksteube 1312 my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
103     mpi_size = mesh_p->Nodes->MPIInfo->size;
104     nDim = mesh_p->Nodes->numDim;
105 dhawcroft 796
106 ksteube 1312 if (! ( (nDim ==2) || (nDim == 3) ) ) {
107     Finley_setError(IO_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
108     return;
109 dhawcroft 793 }
110 ksteube 1312 /*************************************************************************************/
111 dhawcroft 793
112 ksteube 1312 /* open the file and check handle */
113 dhawcroft 793
114 ksteube 1312 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 ksteube 1347 perror(filename_p);
133     sprintf(error_msg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
134 ksteube 1312 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 dhawcroft 793 }
146 ksteube 1312 if (! Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo) ) return;
147     /*************************************************************************************/
148 dhawcroft 793
149 ksteube 1312 /* find the mesh type to be written */
150 dhawcroft 793
151 ksteube 1312 isCellCentered=TMPMEMALLOC(num_data,bool_t);
152     max_len_names=0;
153     if (!Finley_checkPtr(isCellCentered)) {
154 gross 1338 nodetype=FINLEY_UNKNOWN;
155     elementtype=FINLEY_UNKNOWN;
156 ksteube 1312 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 gross 1338 nodetype = (nodetype == FINLEY_UNKNOWN) ? FINLEY_NODES : nodetype;
234 dhawcroft 793 }
235 ksteube 1312 if (Finley_noError()) {
236 dhawcroft 793
237 ksteube 1312 /***************************************/
238 dhawcroft 793
239 ksteube 1312 /* 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 dhawcroft 793 }
393 ksteube 1312 /***************************************/
394 dhawcroft 793
395 ksteube 1312 /***************************************/
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 dhawcroft 793
413 ksteube 1312 /* select number of points and the mesh component */
414 dhawcroft 793
415 ksteube 1312 sprintf(txt_buffer,tags_header,globalNumPoints,numCellFactor*globalNumCells,3);
416 dhawcroft 793
417 ksteube 1312 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 dhawcroft 793
428 ksteube 1312 /* write the nodes */
429    
430     if (mpi_size > 1) {
431 dhawcroft 793
432 ksteube 1312 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 dhawcroft 796
482 ksteube 1312 /* close the Points and open connectivity */
483 dhawcroft 793
484 ksteube 1312 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 dhawcroft 793 }
494    
495 ksteube 1312 /* 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 dhawcroft 793
506 ksteube 1312 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 dhawcroft 796
550 ksteube 1312 }
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 dhawcroft 793 }
560 ksteube 1312 } else {
561     fprintf(fileHandle_p,tags_End_Conn_and_Start_Offset);
562     }
563 dhawcroft 796
564 ksteube 1312 /* 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 dhawcroft 818 }
580 ksteube 1312
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 dhawcroft 793 }
593 ksteube 1312
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 dhawcroft 793 }
615 ksteube 1312 } else {
616     fprintf(fileHandle_p,tags_End_Type_And_Cells);
617 dhawcroft 793 }
618 ksteube 1312 }
619 dhawcroft 818
620 ksteube 1312 /* 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 dhawcroft 793
632 ksteube 1312 switch(getDataPointRank(data_pp[i_data])) {
633 dhawcroft 793 case 0:
634 ksteube 1312 if (! set_scalar) {
635     strcat(txt_buffer," Scalars=\"");
636     strcat(txt_buffer,names_p[i_data]);
637     strcat(txt_buffer,"\"");
638 dhawcroft 793 set_scalar=TRUE;
639     }
640     break;
641     case 1:
642 ksteube 1312 if (! set_vector) {
643     strcat(txt_buffer," Vectors=\"");
644     strcat(txt_buffer,names_p[i_data]);
645     strcat(txt_buffer,"\"");
646 dhawcroft 793 set_vector=TRUE;
647     }
648     break;
649     case 2:
650 ksteube 1312 if (! set_tensor) {
651     strcat(txt_buffer," Tensors=\"");
652     strcat(txt_buffer,names_p[i_data]);
653     strcat(txt_buffer,"\"");
654 dhawcroft 793 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 ksteube 1312 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 dhawcroft 793 }
672 ksteube 1312 } 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 gross 903 } else {
694 ksteube 1312 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 gross 903 }
700 ksteube 1312 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 dhawcroft 793 }
776 ksteube 1312 }
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 dhawcroft 793 }
785 ksteube 1312 } else {
786     fprintf(fileHandle_p,tag_End_CellData);
787 dhawcroft 793 }
788     }
789 ksteube 1312 /* 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 dhawcroft 793
801 ksteube 1312 switch(getDataPointRank(data_pp[i_data])) {
802 dhawcroft 793 case 0:
803 ksteube 1312 if (! set_scalar) {
804     strcat(txt_buffer," Scalars=\"");
805     strcat(txt_buffer,names_p[i_data]);
806     strcat(txt_buffer,"\"");
807 dhawcroft 793 set_scalar=TRUE;
808     }
809     break;
810     case 1:
811 ksteube 1312 if (! set_vector) {
812     strcat(txt_buffer," Vectors=\"");
813     strcat(txt_buffer,names_p[i_data]);
814     strcat(txt_buffer,"\"");
815 dhawcroft 793 set_vector=TRUE;
816     }
817     break;
818     case 2:
819 ksteube 1312 if (! set_tensor) {
820     strcat(txt_buffer," Tensors=\"");
821     strcat(txt_buffer,names_p[i_data]);
822     strcat(txt_buffer,"\"");
823 dhawcroft 793 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 ksteube 1312 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 gross 1062 }
841 ksteube 1312 } 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 dhawcroft 818 }
857 ksteube 1312 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 dhawcroft 818 }
866 ksteube 1312 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 dhawcroft 818 }
872 ksteube 1312 nCompReqd = 9;
873 dhawcroft 793 }
874 ksteube 1312 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 dhawcroft 793 }
941     }
942     }
943 ksteube 1312 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 dhawcroft 793 }
950 gross 1062 } else {
951 ksteube 1312 fprintf(fileHandle_p,tag_End_PointData);
952 gross 1062 }
953 jgs 110 }
954 ksteube 1312 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 dhawcroft 793 }
967 ksteube 1312 } else {
968     fprintf(fileHandle_p,footer);
969     fclose(fileHandle_p);
970     }
971 jgs 153 }
972 gross 1028 TMPMEMFREE(isCellCentered);
973 ksteube 1312 TMPMEMFREE(txt_buffer);
974 jgs 110 return;
975     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26