/[escript]/branches/domexper/dudley/src/Mesh_saveVTK.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3152 - (hide annotations)
Fri Sep 3 05:48:31 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 39341 byte(s)
Collapsed duplicate fields
1 jgs 110
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 caltinay 2141 /***************************************************************************/
15     /* Writes data and mesh in VTK XML format to a VTU file. */
16 jfenwick 3086 /* Nodal data needs to be given on DUDLEY_NODES or DUDLEY_REDUCED_NODES */
17 caltinay 2141 /***************************************************************************/
18 ksteube 1811
19 jgs 110 #include "Mesh.h"
20 gross 1062 #include "Assemble.h"
21 caltinay 2141 #include "vtkCellType.h" /* copied from vtk source directory */
22 gross 1669 #include "paso/PasoUtil.h"
23 jgs 110
24 ksteube 1312 #define INT_FORMAT "%d "
25 caltinay 2744 #define LEN_INT_FORMAT (unsigned int)(9+2)
26 ksteube 1312 #define INT_NEWLINE_FORMAT "%d\n"
27 caltinay 2141 #define SCALAR_FORMAT "%12.6e\n"
28     #define VECTOR_FORMAT "%12.6e %12.6e %12.6e\n"
29     #define TENSOR_FORMAT "%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n"
30 caltinay 2744 /* strlen("-1.234567e+789 ") == 15 */
31     #define LEN_TENSOR_FORMAT (unsigned int)(9*15+2)
32 ksteube 1312 #define NEWLINE "\n"
33 caltinay 2141 #define LEN_TMP_BUFFER LEN_TENSOR_FORMAT+(MAX_numNodes*LEN_INT_FORMAT+1)+2
34 phornby 2172 #define NCOMP_MAX (unsigned int)9
35 caltinay 2141
36     #define __STRCAT(dest, chunk, dest_in_use) \
37     do {\
38     strcpy(&dest[dest_in_use], chunk);\
39     dest_in_use += strlen(chunk);\
40     } while(0)
41    
42     #ifdef PASO_MPI
43     /* writes buffer to file catching the empty buffer case which causes problems
44     * with some MPI versions */
45 gross 2385 #define MPI_WRITE_ORDERED(BUF) \
46 caltinay 2141 do {\
47 gross 2385 int LLEN=0; \
48     LLEN=(int) strlen(BUF); \
49     if (LLEN==0) { strcpy(BUF, ""); LLEN=0; }\
50     MPI_File_write_ordered(mpi_fileHandle_p, BUF, LLEN, MPI_CHAR, &mpi_status);\
51 caltinay 2141 } while(0)
52    
53     /* writes buffer to file on master only */
54     #define MPI_RANK0_WRITE_SHARED(BUF) \
55     do {\
56 gross 2385 int LLEN=0; \
57 caltinay 2141 if (my_mpi_rank == 0) {\
58 gross 2385 LLEN=(int) strlen(BUF); \
59     if (LLEN==0) { strcpy(BUF,""); LLEN=0; }\
60     MPI_File_iwrite_shared(mpi_fileHandle_p, BUF, LLEN, MPI_CHAR, &mpi_req);\
61 caltinay 2141 MPI_Wait(&mpi_req, &mpi_status);\
62     }\
63     } while(0)
64    
65     /* For reference only. Investigation needed as to which values may improve
66     * performance */
67     #if 0
68     void create_MPIInfo(MPI_Info& info)
69     {
70     MPI_Info_create(&info);
71     MPI_Info_set(info, "access_style", "write_once, sequential");
72     MPI_Info_set(info, "collective_buffering", "true");
73     MPI_Info_set(info, "cb_block_size", "131072");
74     MPI_Info_set(info, "cb_buffer_size", "1048567");
75     MPI_Info_set(info, "cb_nodes", "8");
76     MPI_Info_set(info, "striping_factor", "16");
77     MPI_Info_set(info, "striping_unit", "424288");
78 dhawcroft 793 }
79 caltinay 2141 #endif
80 dhawcroft 793
81 caltinay 2141 #else
82    
83 gross 2385 #define MPI_WRITE_ORDERED(A)
84 caltinay 2141 #define MPI_RANK0_WRITE_SHARED(A)
85    
86     #endif /* PASO_MPI */
87    
88    
89     /* Returns one if the node given by coords and idx is within the quadrant
90     * indexed by q and if the element type is Rec9 or Hex27, zero otherwise */
91     int nodeInQuadrant(const double *coords, ElementTypeId type, int idx, int q)
92     {
93     #define INSIDE_1D(_X_,_C_,_R_) ( ABS((_X_)-(_C_)) <= (_R_) )
94     #define INSIDE_2D(_X_,_Y_,_CX_,_CY_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_))
95     #define INSIDE_3D(_X_,_Y_,_Z_,_CX_,_CY_,_CZ_,_R_) ( INSIDE_1D(_X_,_CX_,_R_) && INSIDE_1D(_Y_,_CY_,_R_) && INSIDE_1D(_Z_,_CZ_,_R_) )
96 jfenwick 3126 return 1;
97 caltinay 2141 }
98    
99 jfenwick 3086 void Dudley_Mesh_saveVTK(const char *filename_p,
100     Dudley_Mesh *mesh_p,
101 ksteube 1312 const dim_t num_data,
102 caltinay 2141 char **names_p,
103 gross 2421 escriptDataC **data_pp,
104     const char* metadata,
105     const char*metadata_schema)
106 dhawcroft 793 {
107 caltinay 2141 #ifdef PASO_MPI
108     MPI_File mpi_fileHandle_p;
109     MPI_Status mpi_status;
110     MPI_Request mpi_req;
111     MPI_Info mpi_info = MPI_INFO_NULL;
112     #endif
113     Paso_MPI_rank my_mpi_rank;
114     FILE *fileHandle_p = NULL;
115     char errorMsg[LenErrorMsg_MAX], *txtBuffer;
116     char tmpBuffer[LEN_TMP_BUFFER];
117     size_t txtBufferSize, txtBufferInUse, maxNameLen;
118     double *quadNodes_p = NULL;
119     dim_t dataIdx, nDim;
120     dim_t numCells=0, globalNumCells=0, numVTKNodesPerElement=0;
121     dim_t myNumPoints=0, globalNumPoints=0;
122     dim_t shape, NN=0, numCellFactor=1, myNumCells=0;
123     bool_t *isCellCentered;
124     bool_t writeCellData=FALSE, writePointData=FALSE, hasReducedElements=FALSE;
125     index_t myFirstNode=0, myLastNode=0, *globalNodeIndex=NULL;
126     index_t myFirstCell=0, k;
127     int mpi_size, i, j, l;
128 jfenwick 3086 int cellType=0, nodeType=DUDLEY_NODES, elementType=DUDLEY_UNKNOWN;
129     Dudley_ElementFile *elements = NULL;
130 gross 2748 ElementTypeId typeId = NoRef;
131 dhawcroft 793
132 caltinay 2141 const char *vtkHeader = \
133     "<?xml version=\"1.0\"?>\n" \
134 gross 2421 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"%s%s>\n%s%s" \
135 caltinay 2141 "<UnstructuredGrid>\n" \
136     "<Piece NumberOfPoints=\"%d\" NumberOfCells=\"%d\">\n" \
137     "<Points>\n" \
138     "<DataArray NumberOfComponents=\"%d\" type=\"Float64\" format=\"ascii\">\n";
139     char *vtkFooter = "</Piece>\n</UnstructuredGrid>\n</VTKFile>\n";
140     const char *tag_Float_DataArray="<DataArray Name=\"%s\" type=\"Float64\" NumberOfComponents=\"%d\" format=\"ascii\">\n";
141     char *tags_End_Points_and_Start_Conn = "</DataArray>\n</Points>\n<Cells>\n<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">\n" ;
142     char *tags_End_Conn_and_Start_Offset = "</DataArray>\n<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">\n";
143     char *tags_End_Offset_and_Start_Type = "</DataArray>\n<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n";
144     char *tag_End_DataArray = "</DataArray>\n";
145 dhawcroft 793
146 caltinay 2141 /* if there is no mesh we just return */
147     if (mesh_p==NULL) return;
148 dhawcroft 793
149 caltinay 2141 nDim = mesh_p->Nodes->numDim;
150 dhawcroft 796
151 caltinay 2141 if (nDim != 2 && nDim != 3) {
152 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: spatial dimension 2 or 3 is supported only.");
153 caltinay 2141 return;
154     }
155     my_mpi_rank = mesh_p->Nodes->MPIInfo->rank;
156     mpi_size = mesh_p->Nodes->MPIInfo->size;
157 dhawcroft 793
158 gross 2308 /************************************************************************
159     * open the file and check handle *
160     */
161 caltinay 2141 if (mpi_size > 1) {
162     #ifdef PASO_MPI
163     const int amode = MPI_MODE_CREATE|MPI_MODE_WRONLY|MPI_MODE_UNIQUE_OPEN;
164     int ierr;
165     if (my_mpi_rank == 0 && Paso_fileExists(filename_p)) {
166     remove(filename_p);
167     }
168     ierr = MPI_File_open(mesh_p->Nodes->MPIInfo->comm, (char*)filename_p,
169     amode, mpi_info, &mpi_fileHandle_p);
170     if (ierr != MPI_SUCCESS) {
171     sprintf(errorMsg, "saveVTK: File %s could not be opened for writing in parallel.", filename_p);
172 jfenwick 3086 Dudley_setError(IO_ERROR, errorMsg);
173 caltinay 2141 } else {
174 gross 2310 ierr=MPI_File_set_view(mpi_fileHandle_p,MPI_DISPLACEMENT_CURRENT,
175 caltinay 2141 MPI_CHAR, MPI_CHAR, "native", mpi_info);
176     }
177     #endif /* PASO_MPI */
178     } else {
179 ksteube 1312 fileHandle_p = fopen(filename_p, "w");
180     if (fileHandle_p==NULL) {
181 caltinay 2141 sprintf(errorMsg, "saveVTK: File %s could not be opened for writing.", filename_p);
182 jfenwick 3086 Dudley_setError(IO_ERROR, errorMsg);
183 caltinay 2141 }
184     }
185     if (!Paso_MPIInfo_noError(mesh_p->Nodes->MPIInfo)) return;
186 dhawcroft 793
187 jfenwick 3086 /* General note: From this point if an error occurs Dudley_setError is
188 caltinay 2141 * called and subsequent steps are skipped until the end of this function
189     * where allocated memory is freed and the file is closed. */
190 dhawcroft 793
191 caltinay 2141 /************************************************************************/
192     /* find the mesh type to be written */
193 dhawcroft 793
194 caltinay 2141 isCellCentered = TMPMEMALLOC(num_data, bool_t);
195     maxNameLen = 0;
196 jfenwick 3086 if (!Dudley_checkPtr(isCellCentered)) {
197 caltinay 2141 for (dataIdx=0; dataIdx<num_data; ++dataIdx) {
198     if (! isEmpty(data_pp[dataIdx])) {
199     switch(getFunctionSpaceType(data_pp[dataIdx]) ) {
200 jfenwick 3086 case DUDLEY_NODES:
201     nodeType = (nodeType == DUDLEY_REDUCED_NODES) ? DUDLEY_REDUCED_NODES : DUDLEY_NODES;
202 caltinay 2141 isCellCentered[dataIdx] = FALSE;
203 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
204     elementType = DUDLEY_ELEMENTS;
205 caltinay 2141 } else {
206 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
207 caltinay 2141 }
208     break;
209 jfenwick 3086 case DUDLEY_REDUCED_NODES:
210     nodeType = DUDLEY_REDUCED_NODES;
211 caltinay 2141 isCellCentered[dataIdx] = FALSE;
212 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
213     elementType = DUDLEY_ELEMENTS;
214 caltinay 2141 } else {
215 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
216 caltinay 2141 }
217     break;
218 jfenwick 3086 case DUDLEY_REDUCED_ELEMENTS:
219 caltinay 2141 hasReducedElements = TRUE;
220 jfenwick 3086 case DUDLEY_ELEMENTS:
221 caltinay 2141 isCellCentered[dataIdx] = TRUE;
222 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_ELEMENTS) {
223     elementType = DUDLEY_ELEMENTS;
224 caltinay 2141 } else {
225 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
226 caltinay 2141 }
227     break;
228 jfenwick 3086 case DUDLEY_REDUCED_FACE_ELEMENTS:
229 caltinay 2141 hasReducedElements = TRUE;
230 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
231 caltinay 2141 isCellCentered[dataIdx] = TRUE;
232 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_FACE_ELEMENTS) {
233     elementType = DUDLEY_FACE_ELEMENTS;
234 caltinay 2141 } else {
235 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
236 caltinay 2141 }
237     break;
238 jfenwick 3086 case DUDLEY_POINTS:
239 caltinay 2141 isCellCentered[dataIdx]=TRUE;
240 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN || elementType==DUDLEY_POINTS) {
241     elementType = DUDLEY_POINTS;
242 caltinay 2141 } else {
243 jfenwick 3086 Dudley_setError(TYPE_ERROR, "saveVTK: cannot write given data in single file.");
244 caltinay 2141 }
245     break;
246     default:
247     sprintf(errorMsg, "saveVTK: unknown function space type %d",getFunctionSpaceType(data_pp[dataIdx]));
248 jfenwick 3086 Dudley_setError(TYPE_ERROR, errorMsg);
249 caltinay 2141 }
250     if (isCellCentered[dataIdx]) {
251     writeCellData = TRUE;
252     } else {
253     writePointData = TRUE;
254     }
255     maxNameLen = MAX(maxNameLen, strlen(names_p[dataIdx]));
256     }
257     }
258     }
259 dhawcroft 793
260 caltinay 2141 /************************************************************************/
261     /* select number of points and the mesh component */
262 gross 1741
263 jfenwick 3086 if (Dudley_noError()) {
264     if (nodeType == DUDLEY_REDUCED_NODES) {
265     myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
266     myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
267     globalNumPoints = Dudley_NodeFile_getGlobalNumReducedNodes(mesh_p->Nodes);
268     globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
269 caltinay 2141 } else {
270 jfenwick 3086 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
271     myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
272     globalNumPoints = Dudley_NodeFile_getGlobalNumNodes(mesh_p->Nodes);
273     globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
274 ksteube 1312 }
275 caltinay 2141 myNumPoints = myLastNode - myFirstNode;
276 jfenwick 3086 if (elementType==DUDLEY_UNKNOWN) elementType=DUDLEY_ELEMENTS;
277 caltinay 2141 switch(elementType) {
278 jfenwick 3086 case DUDLEY_ELEMENTS:
279 caltinay 2141 elements = mesh_p->Elements;
280     break;
281 jfenwick 3086 case DUDLEY_FACE_ELEMENTS:
282 caltinay 2141 elements = mesh_p->FaceElements;
283     break;
284 jfenwick 3086 case DUDLEY_POINTS:
285 caltinay 2141 elements = mesh_p->Points;
286     break;
287     }
288     if (elements==NULL) {
289 jfenwick 3086 Dudley_setError(SYSTEM_ERROR, "saveVTK: undefined element file");
290 caltinay 2141 } else {
291 jfenwick 3086 /* map dudley element type to VTK element type */
292 caltinay 2141 numCells = elements->numElements;
293 jfenwick 3086 globalNumCells = Dudley_ElementFile_getGlobalNumElements(elements);
294     myNumCells = Dudley_ElementFile_getMyNumElements(elements);
295     myFirstCell = Dudley_ElementFile_getFirstElement(elements);
296 caltinay 2141 NN = elements->numNodes;
297 caltinay 2814 if (hasReducedElements) {
298 jfenwick 3152 quadNodes_p=elements->referenceElementSet->referenceElementReducedQuadrature->BasisFunctions->QuadNodes;
299 caltinay 2814 } else {
300 jfenwick 3152 quadNodes_p=elements->referenceElementSet->referenceElement->BasisFunctions->QuadNodes;
301 caltinay 2814 }
302 jfenwick 3086 if (nodeType==DUDLEY_REDUCED_NODES) {
303 jfenwick 3152 typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
304 caltinay 2141 } else {
305 gross 2748 typeId = elements->referenceElementSet->referenceElement->Type->TypeId;
306 caltinay 2141 }
307     switch (typeId) {
308     case Point1:
309     case Line2Face:
310     cellType = VTK_VERTEX;
311     numVTKNodesPerElement = 1;
312     break;
313 dhawcroft 793
314 caltinay 2141 case Line2:
315     case Tri3Face:
316     cellType = VTK_LINE;
317     numVTKNodesPerElement = 2;
318     break;
319 dhawcroft 793
320 caltinay 2141 case Tri3:
321     case Tet4Face:
322     cellType = VTK_TRIANGLE;
323     numVTKNodesPerElement = 3;
324     break;
325 dhawcroft 793
326 caltinay 2141 case Tet4:
327     cellType = VTK_TETRA;
328     numVTKNodesPerElement = 4;
329     break;
330 dhawcroft 793
331 caltinay 2141 default:
332 gross 2748 sprintf(errorMsg, "saveVTK: Element type %s is not supported by VTK.", elements->referenceElementSet->referenceElement->Type->Name);
333 jfenwick 3086 Dudley_setError(VALUE_ERROR, errorMsg);
334 caltinay 2141 }
335     }
336     }
337 gross 1741
338 caltinay 2141 /* allocate enough memory for text buffer */
339 gross 1741
340 gross 2421 txtBufferSize = strlen(vtkHeader) + 3*LEN_INT_FORMAT + (30+3*maxNameLen)+strlen(metadata)+strlen(metadata_schema);
341 caltinay 2141 if (mpi_size > 1) {
342 caltinay 2743 txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TMP_BUFFER);
343 caltinay 2141 txtBufferSize = MAX(txtBufferSize, numCellFactor * myNumCells *
344     (LEN_INT_FORMAT * numVTKNodesPerElement + 1));
345     txtBufferSize = MAX(txtBufferSize,
346     numCellFactor * myNumCells * LEN_TENSOR_FORMAT);
347     txtBufferSize = MAX(txtBufferSize, myNumPoints * LEN_TENSOR_FORMAT);
348     }
349     txtBuffer = TMPMEMALLOC(txtBufferSize+1, char);
350 gross 1741
351 caltinay 2141 /* sets error if memory allocation failed */
352 jfenwick 3086 Dudley_checkPtr(txtBuffer);
353 gross 1741
354 caltinay 2141 /************************************************************************/
355     /* write number of points and the mesh component */
356 gross 1741
357 jfenwick 3086 if (Dudley_noError()) {
358 jfenwick 3145 int flag1=0;
359 jfenwick 3086 if (DUDLEY_REDUCED_NODES == nodeType) {
360 jfenwick 3145 flag1=1;
361 jfenwick 3126 } else if (numVTKNodesPerElement != elements->referenceElementSet->referenceElement->Type->numNodes) {
362 jfenwick 3145 flag1=1;
363     }
364 gross 2421 if (strlen(metadata)>0) {
365     if (strlen(metadata_schema)>0) {
366     sprintf(txtBuffer, vtkHeader," ",metadata_schema,metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
367     } else {
368     sprintf(txtBuffer, vtkHeader,"","",metadata,"\n",globalNumPoints, numCellFactor*globalNumCells, 3);
369     }
370     } else {
371     if (strlen(metadata_schema)>0) {
372     sprintf(txtBuffer, vtkHeader," ",metadata_schema,"","",globalNumPoints, numCellFactor*globalNumCells, 3);
373     } else {
374     sprintf(txtBuffer, vtkHeader,"","","","",globalNumPoints, numCellFactor*globalNumCells, 3);
375     }
376     }
377 gross 1741
378 caltinay 2141 if (mpi_size > 1) {
379     /* write the nodes */
380     MPI_RANK0_WRITE_SHARED(txtBuffer);
381     txtBuffer[0] = '\0';
382     txtBufferInUse = 0;
383     if (nDim==2) {
384     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
385     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
386     sprintf(tmpBuffer, VECTOR_FORMAT,
387     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
388     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
389     0.);
390     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
391     }
392     }
393     } else {
394     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
395     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
396     sprintf(tmpBuffer, VECTOR_FORMAT,
397     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
398     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
399     mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
400     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
401     }
402     }
403     } /* nDim */
404 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
405 gross 1741
406 caltinay 2141 /* write the cells */
407     MPI_RANK0_WRITE_SHARED(tags_End_Points_and_Start_Conn);
408     txtBuffer[0] = '\0';
409     txtBufferInUse = 0;
410 jfenwick 3145 if (!flag1) {
411 caltinay 2141 for (i = 0; i < numCells; i++) {
412     if (elements->Owner[i] == my_mpi_rank) {
413     for (j = 0; j < numVTKNodesPerElement; j++) {
414     sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
415     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
416     }
417     __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
418     }
419     }
420     } else {
421     for (i = 0; i < numCells; i++) {
422     if (elements->Owner[i] == my_mpi_rank) {
423     for (l = 0; l < numCellFactor; l++) {
424 jfenwick 3145 const int idx=l*numVTKNodesPerElement;
425 caltinay 2141 for (j = 0; j < numVTKNodesPerElement; j++) {
426 jfenwick 3145 sprintf(tmpBuffer, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx+j, i, NN)]]);
427 caltinay 2141 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
428     }
429     __STRCAT(txtBuffer, NEWLINE, txtBufferInUse);
430     }
431     }
432     }
433     } /* nodeIndex */
434 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
435 gross 1741
436 caltinay 2141 /* write the offsets */
437     MPI_RANK0_WRITE_SHARED(tags_End_Conn_and_Start_Offset);
438     txtBuffer[0] = '\0';
439     txtBufferInUse = 0;
440     for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
441     i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
442     {
443     sprintf(tmpBuffer, INT_NEWLINE_FORMAT, i);
444     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
445     }
446 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
447 gross 1741
448 caltinay 2141 /* write element type */
449     sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
450     MPI_RANK0_WRITE_SHARED(tags_End_Offset_and_Start_Type);
451     txtBuffer[0] = '\0';
452     txtBufferInUse = 0;
453     for (i = numVTKNodesPerElement*(myFirstCell*numCellFactor+1);
454     i <= (myFirstCell+myNumCells)*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement)
455     {
456     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
457     }
458 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
459 caltinay 2141 /* finalize cell information */
460     strcpy(txtBuffer, "</DataArray>\n</Cells>\n");
461     MPI_RANK0_WRITE_SHARED(txtBuffer);
462 gross 1741
463 caltinay 2141 } else { /***** mpi_size == 1 *****/
464 gross 1741
465 caltinay 2141 /* write the nodes */
466     fputs(txtBuffer, fileHandle_p);
467     if (nDim==2) {
468     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
469     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
470     fprintf(fileHandle_p, VECTOR_FORMAT,
471     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
472     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
473     0.);
474     }
475     }
476     } else {
477     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
478     if ( (myFirstNode <= globalNodeIndex[i]) && (globalNodeIndex[i] < myLastNode) ) {
479     fprintf(fileHandle_p, VECTOR_FORMAT,
480     mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)],
481     mesh_p->Nodes->Coordinates[INDEX2(1, i, nDim)],
482     mesh_p->Nodes->Coordinates[INDEX2(2, i, nDim)]);
483     }
484     }
485     } /* nDim */
486 gross 1741
487 caltinay 2141 /* write the cells */
488     fputs(tags_End_Points_and_Start_Conn, fileHandle_p);
489 jfenwick 3145 if (!flag1) {
490 caltinay 2141 for (i = 0; i < numCells; i++) {
491     for (j = 0; j < numVTKNodesPerElement; j++) {
492     fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(j, i, NN)]]);
493     }
494     fprintf(fileHandle_p, NEWLINE);
495     }
496     } else {
497     for (i = 0; i < numCells; i++) {
498     for (l = 0; l < numCellFactor; l++) {
499 jfenwick 3145 const int idx=l*numVTKNodesPerElement;
500 caltinay 2141 for (j = 0; j < numVTKNodesPerElement; j++) {
501 jfenwick 3145 fprintf(fileHandle_p, INT_FORMAT, globalNodeIndex[elements->Nodes[INDEX2(idx+j, i, NN)]]);
502 caltinay 2141 }
503     fprintf(fileHandle_p, NEWLINE);
504     }
505     }
506     } /* nodeIndex */
507 gross 1741
508 caltinay 2141 /* write the offsets */
509     fputs(tags_End_Conn_and_Start_Offset, fileHandle_p);
510     for (i = numVTKNodesPerElement; i <= numCells*numVTKNodesPerElement*numCellFactor; i += numVTKNodesPerElement) {
511     fprintf(fileHandle_p, INT_NEWLINE_FORMAT, i);
512     }
513 gross 1741
514 caltinay 2141 /* write element type */
515     sprintf(tmpBuffer, INT_NEWLINE_FORMAT, cellType);
516     fputs(tags_End_Offset_and_Start_Type, fileHandle_p);
517     for (i = 0; i < numCells*numCellFactor; i++)
518     fputs(tmpBuffer, fileHandle_p);
519     /* finalize cell information */
520     fputs("</DataArray>\n</Cells>\n", fileHandle_p);
521     } /* mpi_size */
522 gross 1741
523 jfenwick 3086 } /* Dudley_noError */
524 dhawcroft 796
525 caltinay 2141 /************************************************************************/
526     /* write cell data */
527 dhawcroft 796
528 jfenwick 3086 if (writeCellData && Dudley_noError()) {
529 caltinay 2141 bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
530     /* mark the active data arrays */
531     strcpy(txtBuffer, "<CellData");
532     for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
533     if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
534     /* rank == 0 <--> scalar data */
535     /* rank == 1 <--> vector data */
536     /* rank == 2 <--> tensor data */
537     switch (getDataPointRank(data_pp[dataIdx])) {
538     case 0:
539     if (!set_scalar) {
540     strcat(txtBuffer, " Scalars=\"");
541     strcat(txtBuffer, names_p[dataIdx]);
542     strcat(txtBuffer, "\"");
543     set_scalar = TRUE;
544     }
545     break;
546     case 1:
547     if (!set_vector) {
548     strcat(txtBuffer, " Vectors=\"");
549     strcat(txtBuffer, names_p[dataIdx]);
550     strcat(txtBuffer, "\"");
551     set_vector = TRUE;
552     }
553     break;
554     case 2:
555     if (!set_tensor) {
556     strcat(txtBuffer, " Tensors=\"");
557     strcat(txtBuffer, names_p[dataIdx]);
558     strcat(txtBuffer, "\"");
559     set_tensor = TRUE;
560     }
561     break;
562     default:
563     sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
564 jfenwick 3086 Dudley_setError(VALUE_ERROR, errorMsg);
565 caltinay 2141 }
566     }
567 jfenwick 3086 if (!Dudley_noError())
568 caltinay 2141 break;
569 dhawcroft 818 }
570 dhawcroft 793 }
571 caltinay 2141 /* only continue if no error occurred */
572 jfenwick 3086 if (writeCellData && Dudley_noError()) {
573 caltinay 2141 strcat(txtBuffer, ">\n");
574     if ( mpi_size > 1) {
575     MPI_RANK0_WRITE_SHARED(txtBuffer);
576     } else {
577     fputs(txtBuffer, fileHandle_p);
578 dhawcroft 793 }
579 dhawcroft 818
580 caltinay 2141 /* write the arrays */
581     for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
582     if (!isEmpty(data_pp[dataIdx]) && isCellCentered[dataIdx]) {
583     dim_t numPointsPerSample=getNumDataPointsPerSample(data_pp[dataIdx]);
584     dim_t rank = getDataPointRank(data_pp[dataIdx]);
585     dim_t nComp = getDataPointSize(data_pp[dataIdx]);
586     dim_t nCompReqd = 1; /* number of components mpi_required */
587     if (rank == 0) {
588     nCompReqd = 1;
589     shape = 0;
590     } else if (rank == 1) {
591     nCompReqd = 3;
592     shape = getDataPointShape(data_pp[dataIdx], 0);
593     if (shape > 3) {
594 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
595 caltinay 2141 }
596     } else {
597     nCompReqd = 9;
598     shape = getDataPointShape(data_pp[dataIdx], 0);
599     if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
600 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
601 caltinay 2141 }
602     }
603     /* bail out if an error occurred */
604 jfenwick 3086 if (!Dudley_noError())
605 caltinay 2141 break;
606 dhawcroft 793
607 caltinay 2141 sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
608     if ( mpi_size > 1) {
609     MPI_RANK0_WRITE_SHARED(txtBuffer);
610     } else {
611     fputs(txtBuffer, fileHandle_p);
612     }
613 gross 1741
614 caltinay 2141 txtBuffer[0] = '\0';
615     txtBufferInUse = 0;
616     for (i=0; i<numCells; i++) {
617     if (elements->Owner[i] == my_mpi_rank) {
618 jfenwick 2770 __const double *values = getSampleDataRO(data_pp[dataIdx], i);
619 caltinay 2141 for (l = 0; l < numCellFactor; l++) {
620     double sampleAvg[NCOMP_MAX];
621     dim_t nCompUsed = MIN(nComp, NCOMP_MAX);
622    
623     /* average over number of points in the sample */
624     if (isExpanded(data_pp[dataIdx])) {
625 gross 2748 dim_t hits=0;
626 caltinay 2141 for (k=0; k<nCompUsed; k++) sampleAvg[k]=0;
627     for (j=0; j<numPointsPerSample; j++) {
628     if (nodeInQuadrant(quadNodes_p, typeId, j, l)) {
629     hits++;
630     for (k=0; k<nCompUsed; k++) {
631     sampleAvg[k] += values[INDEX2(k,j,nComp)];
632     }
633 gross 1741 }
634 caltinay 2141 }
635     for (k=0; k<nCompUsed; k++)
636     sampleAvg[k] /= MAX(hits, 1);
637     } else {
638     for (k=0; k<nCompUsed; k++)
639     sampleAvg[k] = values[k];
640     } /* isExpanded */
641    
642     /* if the number of required components is more than
643     * the number of actual components, pad with zeros
644     */
645     /* probably only need to get shape of first element */
646     if (nCompReqd == 1) {
647     sprintf(tmpBuffer, SCALAR_FORMAT, sampleAvg[0]);
648     } else if (nCompReqd == 3) {
649     if (shape==1) {
650     sprintf(tmpBuffer, VECTOR_FORMAT,
651     sampleAvg[0], 0.f, 0.f);
652     } else if (shape==2) {
653     sprintf(tmpBuffer, VECTOR_FORMAT,
654     sampleAvg[0], sampleAvg[1], 0.f);
655     } else if (shape==3) {
656     sprintf(tmpBuffer, VECTOR_FORMAT,
657     sampleAvg[0],sampleAvg[1],sampleAvg[2]);
658     }
659     } else if (nCompReqd == 9) {
660     if (shape==1) {
661     sprintf(tmpBuffer, TENSOR_FORMAT,
662     sampleAvg[0], 0.f, 0.f,
663     0.f, 0.f, 0.f,
664     0.f, 0.f, 0.f);
665     } else if (shape==2) {
666     sprintf(tmpBuffer, TENSOR_FORMAT,
667     sampleAvg[0], sampleAvg[1], 0.f,
668     sampleAvg[2], sampleAvg[3], 0.f,
669     0.f, 0.f, 0.f);
670     } else if (shape==3) {
671     sprintf(tmpBuffer, TENSOR_FORMAT,
672     sampleAvg[0],sampleAvg[1],sampleAvg[2],
673     sampleAvg[3],sampleAvg[4],sampleAvg[5],
674     sampleAvg[6],sampleAvg[7],sampleAvg[8]);
675     }
676     }
677 gross 1741 if ( mpi_size > 1) {
678 caltinay 2141 __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
679 gross 1741 } else {
680 caltinay 2141 fputs(tmpBuffer, fileHandle_p);
681 gross 1741 }
682 caltinay 2141 } /* for l (numCellFactor) */
683     } /* if I am the owner */
684     } /* for i (numCells) */
685    
686     if ( mpi_size > 1) {
687 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
688 caltinay 2141 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
689     } else {
690     fputs(tag_End_DataArray, fileHandle_p);
691     }
692     } /* !isEmpty && cellCentered */
693     } /* for dataIdx */
694    
695     strcpy(txtBuffer, "</CellData>\n");
696     if ( mpi_size > 1) {
697     MPI_RANK0_WRITE_SHARED(txtBuffer);
698     } else {
699     fputs(txtBuffer, fileHandle_p);
700 dhawcroft 793 }
701 caltinay 2141 } /* if noError && writeCellData */
702 dhawcroft 793
703 caltinay 2141 /************************************************************************/
704     /* write point data */
705    
706 jfenwick 3086 if (writePointData && Dudley_noError()) {
707 caltinay 2141 /* mark the active data arrays */
708     bool_t set_scalar=FALSE, set_vector=FALSE, set_tensor=FALSE;
709     strcpy(txtBuffer, "<PointData");
710     for (dataIdx = 0; dataIdx<num_data; dataIdx++) {
711     if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
712     switch (getDataPointRank(data_pp[dataIdx])) {
713     case 0:
714     if (!set_scalar) {
715     strcat(txtBuffer, " Scalars=\"");
716     strcat(txtBuffer, names_p[dataIdx]);
717     strcat(txtBuffer, "\"");
718     set_scalar = TRUE;
719     }
720     break;
721     case 1:
722     if (!set_vector) {
723     strcat(txtBuffer, " Vectors=\"");
724     strcat(txtBuffer, names_p[dataIdx]);
725     strcat(txtBuffer, "\"");
726     set_vector = TRUE;
727     }
728     break;
729     case 2:
730     if (!set_tensor) {
731     strcat(txtBuffer, " Tensors=\"");
732     strcat(txtBuffer, names_p[dataIdx]);
733     strcat(txtBuffer, "\"");
734     set_tensor = TRUE;
735     }
736     break;
737     default:
738     sprintf(errorMsg, "saveVTK: data %s: VTK supports data with rank <= 2 only.", names_p[dataIdx]);
739 jfenwick 3086 Dudley_setError(VALUE_ERROR, errorMsg);
740 caltinay 2141 }
741 dhawcroft 793 }
742 jfenwick 3086 if (!Dudley_noError())
743 caltinay 2141 break;
744 dhawcroft 793 }
745 caltinay 2141 }
746     /* only continue if no error occurred */
747 jfenwick 3086 if (writePointData && Dudley_noError()) {
748 caltinay 2141 strcat(txtBuffer, ">\n");
749     if ( mpi_size > 1) {
750     MPI_RANK0_WRITE_SHARED(txtBuffer);
751     } else {
752     fputs(txtBuffer, fileHandle_p);
753 gross 1062 }
754 caltinay 2141
755     /* write the arrays */
756     for (dataIdx = 0; dataIdx < num_data; dataIdx++) {
757     if (!isEmpty(data_pp[dataIdx]) && !isCellCentered[dataIdx]) {
758 jfenwick 3086 Dudley_NodeMapping* nodeMapping;
759 caltinay 2141 dim_t rank = getDataPointRank(data_pp[dataIdx]);
760     dim_t nCompReqd = 1; /* number of components mpi_required */
761 jfenwick 3086 if (getFunctionSpaceType(data_pp[dataIdx]) == DUDLEY_REDUCED_NODES) {
762 caltinay 2141 nodeMapping = mesh_p->Nodes->reducedNodesMapping;
763     } else {
764     nodeMapping = mesh_p->Nodes->nodesMapping;
765     }
766     if (rank == 0) {
767     nCompReqd = 1;
768     shape = 0;
769     } else if (rank == 1) {
770     nCompReqd = 3;
771     shape = getDataPointShape(data_pp[dataIdx], 0);
772     if (shape > 3) {
773 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 1 objects must have 3 components at most.");
774 caltinay 2141 }
775     } else {
776     nCompReqd = 9;
777     shape=getDataPointShape(data_pp[dataIdx], 0);
778     if (shape > 3 || shape != getDataPointShape(data_pp[dataIdx], 1)) {
779 jfenwick 3086 Dudley_setError(VALUE_ERROR, "saveVTK: rank 2 objects of shape 2x2 or 3x3 supported only.");
780 caltinay 2141 }
781     }
782     /* bail out if an error occurred */
783 jfenwick 3086 if (!Dudley_noError())
784 caltinay 2141 break;
785    
786     sprintf(txtBuffer, tag_Float_DataArray, names_p[dataIdx], nCompReqd);
787     if ( mpi_size > 1) {
788     MPI_RANK0_WRITE_SHARED(txtBuffer);
789     } else {
790     fputs(txtBuffer, fileHandle_p);
791     }
792    
793     txtBuffer[0] = '\0';
794     txtBufferInUse = 0;
795     for (i=0; i<mesh_p->Nodes->numNodes; i++) {
796     k = globalNodeIndex[i];
797     if ( (myFirstNode <= k) && (k < myLastNode) ) {
798 jfenwick 2770 __const double *values = getSampleDataRO(data_pp[dataIdx], nodeMapping->target[i]);
799 caltinay 2141 /* if the number of mpi_required components is more than
800     * the number of actual components, pad with zeros.
801     * Probably only need to get shape of first element */
802     if (nCompReqd == 1) {
803     sprintf(tmpBuffer, SCALAR_FORMAT, values[0]);
804     } else if (nCompReqd == 3) {
805     if (shape==1) {
806     sprintf(tmpBuffer, VECTOR_FORMAT,
807     values[0], 0.f, 0.f);
808     } else if (shape==2) {
809     sprintf(tmpBuffer, VECTOR_FORMAT,
810     values[0], values[1], 0.f);
811     } else if (shape==3) {
812     sprintf(tmpBuffer, VECTOR_FORMAT,
813     values[0], values[1], values[2]);
814     }
815     } else if (nCompReqd == 9) {
816     if (shape==1) {
817     sprintf(tmpBuffer, TENSOR_FORMAT,
818     values[0], 0.f, 0.f,
819     0.f, 0.f, 0.f,
820     0.f, 0.f, 0.f);
821     } else if (shape==2) {
822     sprintf(tmpBuffer, TENSOR_FORMAT,
823     values[0], values[1], 0.f,
824     values[2], values[3], 0.f,
825     0.f, 0.f, 0.f);
826     } else if (shape==3) {
827     sprintf(tmpBuffer, TENSOR_FORMAT,
828     values[0], values[1], values[2],
829     values[3], values[4], values[5],
830     values[6], values[7], values[8]);
831     }
832     }
833     if ( mpi_size > 1) {
834     __STRCAT(txtBuffer, tmpBuffer, txtBufferInUse);
835     } else {
836     fputs(tmpBuffer, fileHandle_p);
837     }
838     } /* if this is my node */
839     } /* for i (numNodes) */
840    
841     if ( mpi_size > 1) {
842 gross 2385 MPI_WRITE_ORDERED(txtBuffer);
843 caltinay 2141 MPI_RANK0_WRITE_SHARED(tag_End_DataArray);
844     } else {
845     fputs(tag_End_DataArray, fileHandle_p);
846     }
847     } /* !isEmpty && !isCellCentered */
848     } /* for dataIdx */
849    
850     strcpy(txtBuffer, "</PointData>\n");
851 ksteube 1312 if ( mpi_size > 1) {
852 caltinay 2141 MPI_RANK0_WRITE_SHARED(txtBuffer);
853 gross 1062 } else {
854 caltinay 2141 fputs(txtBuffer, fileHandle_p);
855 gross 1062 }
856 caltinay 2141 } /* if noError && writePointData */
857    
858     /* Final write to VTK file */
859 jfenwick 3086 if (Dudley_noError()) {
860 caltinay 2141 if (mpi_size > 1) {
861     MPI_RANK0_WRITE_SHARED(vtkFooter);
862     } else {
863     fputs(vtkFooter, fileHandle_p);
864 dhawcroft 793 }
865 caltinay 2141 }
866    
867     if ( mpi_size > 1) {
868 caltinay 2094 #ifdef PASO_MPI
869 caltinay 2141 MPI_File_close(&mpi_fileHandle_p);
870 caltinay 2743 MPI_Barrier(mesh_p->Nodes->MPIInfo->comm);
871 caltinay 2094 #endif
872 caltinay 2141 } else {
873     fclose(fileHandle_p);
874     }
875     TMPMEMFREE(isCellCentered);
876     TMPMEMFREE(txtBuffer);
877 jgs 110 }
878 caltinay 2141

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26