/[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 3220 - (hide annotations)
Wed Sep 29 00:33:16 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 39647 byte(s)
Removing references to Quadrature.? and ShapeFns

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26