Actual source code: grvtk.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/daimpl.h>
2: #include <../src/sys/viewer/impls/vtk/vtkvimpl.h>
4: #if defined(PETSC_HAVE_STDINT_H) /* The VTK format requires a 32-bit integer */
5: typedef int32_t PetscVTKInt;
6: #else
7: typedef int PetscVTKInt;
8: #endif
9: #define PETSC_VTK_INT_MAX 2147483647
10: #define PETSC_VTK_INT_MIN -2147483647
11: #if defined(PETSC_USE_64BIT_INDICES)
12: # define PetscVTKIntCheck(a) if ((a) > PETSC_VTK_INT_MAX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Array too long for 32-bit VTK binary format")
13: # define PetscVTKIntCast(a) (PetscVTKInt)(a);PetscVTKIntCheck(a)
14: #else
15: # define PetscVTKIntCheck(a)
16: # define PetscVTKIntCast(a) a
17: #endif
21: /* Write binary data preceded by 32-bit int length (in bytes), does not do byte swapping. */
22: static PetscErrorCode PetscFWrite_VTK(MPI_Comm comm,FILE *fp,void *data,PetscInt n,PetscDataType dtype)
23: {
25: PetscMPIInt rank;
28: if (n < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Trying to write a negative amount of data %D",n);
29: if (!n) return(0);
30: MPI_Comm_rank(comm,&rank);
31: if (!rank) {
32: size_t count;
33: PetscInt size;
34: PetscVTKInt bytes;
35: switch (dtype) {
36: case PETSC_DOUBLE:
37: size = sizeof(double);
38: break;
39: case PETSC_FLOAT:
40: size = sizeof(float);
41: break;
42: case PETSC_INT:
43: size = sizeof(PetscInt);
44: break;
45: case PETSC_CHAR:
46: size = sizeof(char);
47: break;
48: default: SETERRQ(comm,PETSC_ERR_SUP,"Data type not supported");
49: }
50: bytes = PetscVTKIntCast(size*n);
52: count = fwrite(&bytes,sizeof(int),1,fp);
53: if (count != 1) {
54: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_WRITE,"Error writing byte count");
55: }
56: count = fwrite(data,size,(size_t)n,fp);
57: if ((PetscInt)count != n) {
58: SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_FILE_WRITE,"Wrote %D/%D array members of size %D",(PetscInt)count,n,(PetscInt)size);
59: }
60: }
61: return(0);
62: }
66: static PetscErrorCode DMDAVTKWriteAll_VTS(DM da,PetscViewer viewer)
67: {
68: #if defined(PETSC_USE_REAL_SINGLE)
69: const char precision[] = "Float32";
70: #elif defined(PETSC_USE_REAL_DOUBLE)
71: const char precision[] = "Float64";
72: #else
73: const char precision[] = "UnknownPrecision";
74: #endif
75: MPI_Comm comm = ((PetscObject)da)->comm;
76: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
77: PetscViewerVTKObjectLink link;
78: FILE *fp;
79: PetscMPIInt rank,size,tag;
80: DMDALocalInfo info;
81: PetscInt dim,mx,my,mz,bs,boffset,maxnnodes,i,j,k,f,r;
82: PetscInt rloc[6],(*grloc)[6] = PETSC_NULL;
83: PetscScalar *array,*array2;
84: PetscReal gmin[3],gmax[3];
85: PetscErrorCode ierr;
88: #if defined(PETSC_USE_COMPLEX)
89: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Complex values not supported");
90: #endif
92: MPI_Comm_size(comm,&size);
93: MPI_Comm_rank(comm,&rank);
94: DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);
95: DMDAGetLocalInfo(da,&info);
96: DMDAGetBoundingBox(da,gmin,gmax);
98: PetscFOpen(comm,vtk->filename,"wb",&fp);
99: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
100: #ifdef PETSC_WORDS_BIGENDIAN
101: PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
102: #else
103: PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
104: #endif
105: PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);
107: if (!rank) {PetscMalloc(size*6*sizeof(PetscInt),&grloc);}
108: rloc[0] = info.xs;
109: rloc[1] = info.xm;
110: rloc[2] = info.ys;
111: rloc[3] = info.ym;
112: rloc[4] = info.zs;
113: rloc[5] = info.zm;
114: MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);
116: /* Write XML header */
117: maxnnodes = 0; /* Used for the temporary array size on rank 0 */
118: boffset = 0; /* Offset into binary file */
119: for (r=0; r<size; r++) {
120: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
121: if (!rank) {
122: xs = grloc[r][0];
123: xm = grloc[r][1];
124: ys = grloc[r][2];
125: ym = grloc[r][3];
126: zs = grloc[r][4];
127: zm = grloc[r][5];
128: nnodes = xm*ym*zm;
129: }
130: maxnnodes = PetscMax(maxnnodes,nnodes);
131: #if 0
132: switch (dim) {
133: case 1:
134: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);
135: break;
136: case 2:
137: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm,ys+ym-1,xs,xs+xm-1,0,0);
138: break;
139: case 3:
140: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
141: break;
142: default: SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"No support for dimension %D",dim);
143: }
144: #endif
145: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
146: PetscFPrintf(comm,fp," <Points>\n");
147: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
148: boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int);
149: PetscFPrintf(comm,fp," </Points>\n");
151: PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");
152: for (link=vtk->link; link; link=link->next) {
153: Vec X = (Vec)link->vec;
154: const char *vecname = "";
155: if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
156: PetscObjectGetName((PetscObject)X,&vecname);
157: }
158: for (i=0; i<bs; i++) {
159: char buf[256];
160: const char *fieldname;
161: DMDAGetFieldName(da,i,&fieldname);
162: if (!fieldname) {
163: PetscSNPrintf(buf,sizeof buf,"Unnamed%D",i);
164: fieldname = buf;
165: }
166: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
167: boffset += nnodes*sizeof(PetscScalar) + sizeof(int);
168: }
169: }
170: PetscFPrintf(comm,fp," </PointData>\n");
171: PetscFPrintf(comm,fp," </Piece>\n");
172: }
173: PetscFPrintf(comm,fp," </StructuredGrid>\n");
174: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
175: PetscFPrintf(comm,fp,"_");
177: /* Now write the arrays. */
178: tag = ((PetscObject)viewer)->tag;
179: PetscMalloc2(maxnnodes*PetscMax(3,bs),PetscScalar,&array,maxnnodes*3,PetscScalar,&array2);
180: for (r=0; r<size; r++) {
181: MPI_Status status;
182: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
183: if (!rank) {
184: xs = grloc[r][0];
185: xm = grloc[r][1];
186: ys = grloc[r][2];
187: ym = grloc[r][3];
188: zs = grloc[r][4];
189: zm = grloc[r][5];
190: nnodes = xm*ym*zm;
191: } else if (r == rank) {
192: nnodes = info.xm*info.ym*info.zm;
193: }
195: { /* Write the coordinates */
196: Vec Coords;
197: DMDAGetCoordinates(da,&Coords);
198: if (Coords) {
199: const PetscScalar *coords;
200: VecGetArrayRead(Coords,&coords);
201: if (!rank) {
202: if (r) {
203: PetscMPIInt nn;
204: MPI_Recv(array,nnodes*dim,MPIU_SCALAR,r,tag,comm,&status);
205: MPI_Get_count(&status,MPIU_SCALAR,&nn);
206: if (nn != nnodes*dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
207: } else {
208: PetscMemcpy(array,coords,nnodes*dim*sizeof(PetscScalar));
209: }
210: /* Transpose coordinates to VTK (C-style) ordering */
211: for (k=0; k<zm; k++) {
212: for (j=0; j<ym; j++) {
213: for (i=0; i<xm; i++) {
214: PetscInt Iloc = i+xm*(j+ym*k);
215: array2[Iloc*3+0] = array[Iloc*dim + 0];
216: array2[Iloc*3+1] = dim > 1 ? array[Iloc*dim + 1] : 0;
217: array2[Iloc*3+2] = dim > 2 ? array[Iloc*dim + 2] : 0;
218: }
219: }
220: }
221: } else if (r == rank) {
222: MPI_Send((void*)coords,nnodes*dim,MPIU_SCALAR,0,tag,comm);
223: }
224: VecRestoreArrayRead(Coords,&coords);
225: } else { /* Fabricate some coordinates using grid index */
226: for (k=0; k<zm; k++) {
227: for (j=0; j<ym; j++) {
228: for (i=0; i<xm; i++) {
229: PetscInt Iloc = i+xm*(j+ym*k);
230: array2[Iloc*3+0] = xs+i;
231: array2[Iloc*3+1] = ys+j;
232: array2[Iloc*3+2] = zs+k;
233: }
234: }
235: }
236: }
237: PetscFWrite_VTK(comm,fp,array2,nnodes*3,PETSC_SCALAR);
238: }
240: /* Write each of the objects queued up for this file */
241: for (link=vtk->link; link; link=link->next) {
242: Vec X = (Vec)link->vec;
243: const PetscScalar *x;
245: VecGetArrayRead(X,&x);
246: if (!rank) {
247: if (r) {
248: PetscMPIInt nn;
249: MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);
250: MPI_Get_count(&status,MPIU_SCALAR,&nn);
251: if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r);
252: } else {
253: PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));
254: }
255: for (f=0; f<bs; f++) {
256: /* Extract and transpose the f'th field */
257: for (k=0; k<zm; k++) {
258: for (j=0; j<ym; j++) {
259: for (i=0; i<xm; i++) {
260: PetscInt Iloc = i+xm*(j+ym*k);
261: array2[Iloc] = array[Iloc*bs + f];
262: }
263: }
264: }
265: PetscFWrite_VTK(comm,fp,array2,nnodes,PETSC_SCALAR);
266: }
267: } else if (r == rank) {
268: MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);
269: }
270: VecRestoreArrayRead(X,&x);
271: }
272: }
273: PetscFree2(array,array2);
274: PetscFree(grloc);
276: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
277: PetscFPrintf(comm,fp,"</VTKFile>\n");
278: PetscFClose(comm,fp);
279: return(0);
280: }
285: /*@C
286: DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
288: Collective
290: Input Arguments:
291: odm - DM specifying the grid layout, passed as a PetscObject
292: viewer - viewer of type VTK
294: Level: developer
296: Note:
297: This function is a callback used by the VTK viewer to actually write the file.
298: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
299: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
301: .seealso: PETSCVIEWERVTK
302: @*/
303: PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer)
304: {
305: DM dm = (DM)odm;
306: PetscBool isvtk;
307: PetscErrorCode ierr;
312: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);
313: if (!isvtk) SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name);
314: switch (viewer->format) {
315: case PETSC_VIEWER_VTK_VTS:
316: DMDAVTKWriteAll_VTS(dm,viewer);
317: break;
318: default: SETERRQ1(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]);
319: }
320: return(0);
321: }