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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26