Actual source code: grvtk.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/dmdaimpl.h>
2: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
6: static PetscErrorCode DMDAVTKWriteAll_VTS(DM da,PetscViewer viewer)
7: {
8: #if defined(PETSC_USE_REAL_SINGLE)
9: const char precision[] = "Float32";
10: #elif defined(PETSC_USE_REAL_DOUBLE)
11: const char precision[] = "Float64";
12: #else
13: const char precision[] = "UnknownPrecision";
14: #endif
15: MPI_Comm comm;
16: Vec Coords;
17: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
18: PetscViewerVTKObjectLink link;
19: FILE *fp;
20: PetscMPIInt rank,size,tag;
21: DMDALocalInfo info;
22: PetscInt dim,mx,my,mz,cdim,bs,boffset,maxnnodes,i,j,k,f,r;
23: PetscInt rloc[6],(*grloc)[6] = NULL;
24: PetscScalar *array,*array2;
25: PetscReal gmin[3],gmax[3];
26: PetscErrorCode ierr;
29: PetscObjectGetComm((PetscObject)da,&comm);
30: #if defined(PETSC_USE_COMPLEX)
31: SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported");
32: #endif
33: MPI_Comm_size(comm,&size);
34: MPI_Comm_rank(comm,&rank);
35: DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);
36: DMDAGetLocalInfo(da,&info);
37: DMDAGetBoundingBox(da,gmin,gmax);
38: DMGetCoordinates(da,&Coords);
39: if (Coords) {
40: PetscInt csize;
41: VecGetSize(Coords,&csize);
42: if (csize % (mx*my*mz)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate vector size mismatch");
43: cdim = csize/(mx*my*mz);
44: if (cdim < dim || cdim > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate vector size mismatch");
45: } else {
46: cdim = dim;
47: }
49: PetscFOpen(comm,vtk->filename,"wb",&fp);
50: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
51: #if defined(PETSC_WORDS_BIGENDIAN)
52: PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
53: #else
54: PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
55: #endif
56: PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);
58: if (!rank) {PetscMalloc1(size*6,&grloc);}
59: rloc[0] = info.xs;
60: rloc[1] = info.xm;
61: rloc[2] = info.ys;
62: rloc[3] = info.ym;
63: rloc[4] = info.zs;
64: rloc[5] = info.zm;
65: MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);
67: /* Write XML header */
68: maxnnodes = 0; /* Used for the temporary array size on rank 0 */
69: boffset = 0; /* Offset into binary file */
70: for (r=0; r<size; r++) {
71: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
72: if (!rank) {
73: xs = grloc[r][0];
74: xm = grloc[r][1];
75: ys = grloc[r][2];
76: ym = grloc[r][3];
77: zs = grloc[r][4];
78: zm = grloc[r][5];
79: nnodes = xm*ym*zm;
80: }
81: maxnnodes = PetscMax(maxnnodes,nnodes);
82: #if 0
83: switch (dim) {
84: case 1:
85: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);
86: break;
87: case 2:
88: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm,ys+ym-1,xs,xs+xm-1,0,0);
89: break;
90: case 3:
91: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
92: break;
93: default: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"No support for dimension %D",dim);
94: }
95: #endif
96: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
97: PetscFPrintf(comm,fp," <Points>\n");
98: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
99: boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int);
100: PetscFPrintf(comm,fp," </Points>\n");
102: PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");
103: for (link=vtk->link; link; link=link->next) {
104: Vec X = (Vec)link->vec;
105: const char *vecname = "";
106: 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. */
107: PetscObjectGetName((PetscObject)X,&vecname);
108: }
109: for (i=0; i<bs; i++) {
110: char buf[256];
111: const char *fieldname;
112: DMDAGetFieldName(da,i,&fieldname);
113: if (!fieldname) {
114: PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);
115: fieldname = buf;
116: }
117: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
118: boffset += nnodes*sizeof(PetscScalar) + sizeof(int);
119: }
120: }
121: PetscFPrintf(comm,fp," </PointData>\n");
122: PetscFPrintf(comm,fp," </Piece>\n");
123: }
124: PetscFPrintf(comm,fp," </StructuredGrid>\n");
125: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
126: PetscFPrintf(comm,fp,"_");
128: /* Now write the arrays. */
129: tag = ((PetscObject)viewer)->tag;
130: PetscMalloc2(maxnnodes*PetscMax(3,bs),&array,maxnnodes*3,&array2);
131: for (r=0; r<size; r++) {
132: MPI_Status status;
133: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
134: if (!rank) {
135: xs = grloc[r][0];
136: xm = grloc[r][1];
137: ys = grloc[r][2];
138: ym = grloc[r][3];
139: zs = grloc[r][4];
140: zm = grloc[r][5];
141: nnodes = xm*ym*zm;
142: } else if (r == rank) {
143: nnodes = info.xm*info.ym*info.zm;
144: }
146: /* Write the coordinates */
147: if (Coords) {
148: const PetscScalar *coords;
149: VecGetArrayRead(Coords,&coords);
150: if (!rank) {
151: if (r) {
152: PetscMPIInt nn;
153: MPI_Recv(array,nnodes*cdim,MPIU_SCALAR,r,tag,comm,&status);
154: MPI_Get_count(&status,MPIU_SCALAR,&nn);
155: if (nn != nnodes*cdim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
156: } else {
157: PetscMemcpy(array,coords,nnodes*cdim*sizeof(PetscScalar));
158: }
159: /* Transpose coordinates to VTK (C-style) ordering */
160: for (k=0; k<zm; k++) {
161: for (j=0; j<ym; j++) {
162: for (i=0; i<xm; i++) {
163: PetscInt Iloc = i+xm*(j+ym*k);
164: array2[Iloc*3+0] = array[Iloc*cdim + 0];
165: array2[Iloc*3+1] = cdim > 1 ? array[Iloc*cdim + 1] : 0.0;
166: array2[Iloc*3+2] = cdim > 2 ? array[Iloc*cdim + 2] : 0.0;
167: }
168: }
169: }
170: } else if (r == rank) {
171: MPI_Send((void*)coords,nnodes*cdim,MPIU_SCALAR,0,tag,comm);
172: }
173: VecRestoreArrayRead(Coords,&coords);
174: } else { /* Fabricate some coordinates using grid index */
175: for (k=0; k<zm; k++) {
176: for (j=0; j<ym; j++) {
177: for (i=0; i<xm; i++) {
178: PetscInt Iloc = i+xm*(j+ym*k);
179: array2[Iloc*3+0] = xs+i;
180: array2[Iloc*3+1] = ys+j;
181: array2[Iloc*3+2] = zs+k;
182: }
183: }
184: }
185: }
186: PetscViewerVTKFWrite(viewer,fp,array2,nnodes*3,PETSC_SCALAR);
188: /* Write each of the objects queued up for this file */
189: for (link=vtk->link; link; link=link->next) {
190: Vec X = (Vec)link->vec;
191: const PetscScalar *x;
193: VecGetArrayRead(X,&x);
194: if (!rank) {
195: if (r) {
196: PetscMPIInt nn;
197: MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);
198: MPI_Get_count(&status,MPIU_SCALAR,&nn);
199: if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r);
200: } else {
201: PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));
202: }
203: for (f=0; f<bs; f++) {
204: /* Extract and transpose the f'th field */
205: for (k=0; k<zm; k++) {
206: for (j=0; j<ym; j++) {
207: for (i=0; i<xm; i++) {
208: PetscInt Iloc = i+xm*(j+ym*k);
209: array2[Iloc] = array[Iloc*bs + f];
210: }
211: }
212: }
213: PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);
214: }
215: } else if (r == rank) {
216: MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);
217: }
218: VecRestoreArrayRead(X,&x);
219: }
220: }
221: PetscFree2(array,array2);
222: PetscFree(grloc);
224: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
225: PetscFPrintf(comm,fp,"</VTKFile>\n");
226: PetscFClose(comm,fp);
227: return(0);
228: }
233: static PetscErrorCode DMDAVTKWriteAll_VTR(DM da,PetscViewer viewer)
234: {
235: #if defined(PETSC_USE_REAL_SINGLE)
236: const char precision[] = "Float32";
237: #elif defined(PETSC_USE_REAL_DOUBLE)
238: const char precision[] = "Float64";
239: #else
240: const char precision[] = "UnknownPrecision";
241: #endif
242: MPI_Comm comm;
243: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
244: PetscViewerVTKObjectLink link;
245: FILE *fp;
246: PetscMPIInt rank,size,tag;
247: DMDALocalInfo info;
248: PetscInt dim,mx,my,mz,bs,boffset,maxnnodes,i,j,k,f,r;
249: PetscInt rloc[6],(*grloc)[6] = NULL;
250: PetscScalar *array,*array2;
251: PetscReal gmin[3],gmax[3];
252: PetscErrorCode ierr;
255: PetscObjectGetComm((PetscObject)da,&comm);
256: #if defined(PETSC_USE_COMPLEX)
257: SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported");
258: #endif
259: MPI_Comm_size(comm,&size);
260: MPI_Comm_rank(comm,&rank);
261: DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);
262: DMDAGetLocalInfo(da,&info);
263: DMDAGetBoundingBox(da,gmin,gmax);
264: PetscFOpen(comm,vtk->filename,"wb",&fp);
265: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
266: #if defined(PETSC_WORDS_BIGENDIAN)
267: PetscFPrintf(comm,fp,"<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
268: #else
269: PetscFPrintf(comm,fp,"<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
270: #endif
271: PetscFPrintf(comm,fp," <RectilinearGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);
273: if (!rank) {PetscMalloc1(size*6,&grloc);}
274: rloc[0] = info.xs;
275: rloc[1] = info.xm;
276: rloc[2] = info.ys;
277: rloc[3] = info.ym;
278: rloc[4] = info.zs;
279: rloc[5] = info.zm;
280: MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);
282: /* Write XML header */
283: maxnnodes = 0; /* Used for the temporary array size on rank 0 */
284: boffset = 0; /* Offset into binary file */
285: for (r=0; r<size; r++) {
286: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
287: if (!rank) {
288: xs = grloc[r][0];
289: xm = grloc[r][1];
290: ys = grloc[r][2];
291: ym = grloc[r][3];
292: zs = grloc[r][4];
293: zm = grloc[r][5];
294: nnodes = xm*ym*zm;
295: }
296: maxnnodes = PetscMax(maxnnodes,nnodes);
297: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
298: PetscFPrintf(comm,fp," <Coordinates>\n");
299: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Xcoord\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
300: boffset += xm*sizeof(PetscScalar) + sizeof(int);
301: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Ycoord\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
302: boffset += ym*sizeof(PetscScalar) + sizeof(int);
303: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Zcoord\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
304: boffset += zm*sizeof(PetscScalar) + sizeof(int);
305: PetscFPrintf(comm,fp," </Coordinates>\n");
306: PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");
307: for (link=vtk->link; link; link=link->next) {
308: Vec X = (Vec)link->vec;
309: const char *vecname = "";
310: 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. */
311: PetscObjectGetName((PetscObject)X,&vecname);
312: }
313: for (i=0; i<bs; i++) {
314: char buf[256];
315: const char *fieldname;
316: DMDAGetFieldName(da,i,&fieldname);
317: if (!fieldname) {
318: PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);
319: fieldname = buf;
320: }
321: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
322: boffset += nnodes*sizeof(PetscScalar) + sizeof(int);
323: }
324: }
325: PetscFPrintf(comm,fp," </PointData>\n");
326: PetscFPrintf(comm,fp," </Piece>\n");
327: }
328: PetscFPrintf(comm,fp," </RectilinearGrid>\n");
329: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
330: PetscFPrintf(comm,fp,"_");
332: /* Now write the arrays. */
333: tag = ((PetscObject)viewer)->tag;
334: PetscMalloc2(PetscMax(maxnnodes*bs,info.xm+info.ym+info.zm),&array,PetscMax(maxnnodes*3,info.xm+info.ym+info.zm),&array2);
335: for (r=0; r<size; r++) {
336: MPI_Status status;
337: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
338: if (!rank) {
339: xs = grloc[r][0];
340: xm = grloc[r][1];
341: ys = grloc[r][2];
342: ym = grloc[r][3];
343: zs = grloc[r][4];
344: zm = grloc[r][5];
345: nnodes = xm*ym*zm;
346: } else if (r == rank) {
347: nnodes = info.xm*info.ym*info.zm;
348: }
349: { /* Write the coordinates */
350: Vec Coords;
351: DMGetCoordinates(da,&Coords);
352: if (Coords) {
353: const PetscScalar *coords;
354: VecGetArrayRead(Coords,&coords);
355: if (!rank) {
356: if (r) {
357: PetscMPIInt nn;
358: MPI_Recv(array,xm+ym+zm,MPIU_SCALAR,r,tag,comm,&status);
359: MPI_Get_count(&status,MPIU_SCALAR,&nn);
360: if (nn != xm+ym+zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
361: } else {
362: /* x coordinates */
363: for (j=0, k=0, i=0; i<xm; i++) {
364: PetscInt Iloc = i+xm*(j+ym*k);
365: array[i] = coords[Iloc*dim + 0];
366: }
367: /* y coordinates */
368: for (i=0, k=0, j=0; j<ym; j++) {
369: PetscInt Iloc = i+xm*(j+ym*k);
370: array[j+xm] = dim > 1 ? coords[Iloc*dim + 1] : 0;
371: }
372: /* z coordinates */
373: for (i=0, j=0, k=0; k<zm; k++) {
374: PetscInt Iloc = i+xm*(j+ym*k);
375: array[k+xm+ym] = dim > 2 ? coords[Iloc*dim + 2] : 0;
376: }
377: }
378: } else if (r == rank) {
379: xm = info.xm;
380: ym = info.ym;
381: zm = info.zm;
382: for (j=0, k=0, i=0; i<xm; i++) {
383: PetscInt Iloc = i+xm*(j+ym*k);
384: array2[i] = coords[Iloc*dim + 0];
385: }
386: for (i=0, k=0, j=0; j<ym; j++) {
387: PetscInt Iloc = i+xm*(j+ym*k);
388: array2[j+xm] = dim > 1 ? coords[Iloc*dim + 1] : 0;
389: }
390: for (i=0, j=0, k=0; k<zm; k++) {
391: PetscInt Iloc = i+xm*(j+ym*k);
392: array2[k+xm+ym] = dim > 2 ? coords[Iloc*dim + 2] : 0;
393: }
394: MPI_Send((void*)array2,xm+ym+zm,MPIU_SCALAR,0,tag,comm);
395: }
396: VecRestoreArrayRead(Coords,&coords);
397: } else { /* Fabricate some coordinates using grid index */
398: for (j=0, k=0, i=0; i<xm; i++) {
399: array[i] = xs+i;
400: }
401: for (i=0, k=0, j=0; j<ym; j++) {
402: array[j+xm] = ys+j;
403: }
404: for (i=0, j=0, k=0; k<zm; k++) {
405: array[k+xm+ym] = zs+k;
406: }
407: }
408: if (!rank) {
409: PetscViewerVTKFWrite(viewer,fp,&(array[0]),xm,PETSC_SCALAR);
410: PetscViewerVTKFWrite(viewer,fp,&(array[xm]),ym,PETSC_SCALAR);
411: PetscViewerVTKFWrite(viewer,fp,&(array[xm+ym]),zm,PETSC_SCALAR);
412: }
413: }
415: /* Write each of the objects queued up for this file */
416: for (link=vtk->link; link; link=link->next) {
417: Vec X = (Vec)link->vec;
418: const PetscScalar *x;
420: VecGetArrayRead(X,&x);
421: if (!rank) {
422: if (r) {
423: PetscMPIInt nn;
424: MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);
425: MPI_Get_count(&status,MPIU_SCALAR,&nn);
426: if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r);
427: } else {
428: PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));
429: }
430: for (f=0; f<bs; f++) {
431: /* Extract and transpose the f'th field */
432: for (k=0; k<zm; k++) {
433: for (j=0; j<ym; j++) {
434: for (i=0; i<xm; i++) {
435: PetscInt Iloc = i+xm*(j+ym*k);
436: array2[Iloc] = array[Iloc*bs + f];
437: }
438: }
439: }
440: PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);
441: }
442: } else if (r == rank) {
443: MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);
444: }
445: VecRestoreArrayRead(X,&x);
446: }
447: }
448: PetscFree2(array,array2);
449: PetscFree(grloc);
451: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
452: PetscFPrintf(comm,fp,"</VTKFile>\n");
453: PetscFClose(comm,fp);
454: return(0);
455: }
459: /*@C
460: DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
462: Collective
464: Input Arguments:
465: odm - DM specifying the grid layout, passed as a PetscObject
466: viewer - viewer of type VTK
468: Level: developer
470: Note:
471: This function is a callback used by the VTK viewer to actually write the file.
472: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
473: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
475: .seealso: PETSCVIEWERVTK
476: @*/
477: PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer)
478: {
479: DM dm = (DM)odm;
480: PetscBool isvtk;
486: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);
487: if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name);
488: switch (viewer->format) {
489: case PETSC_VIEWER_VTK_VTS:
490: DMDAVTKWriteAll_VTS(dm,viewer);
491: break;
492: case PETSC_VIEWER_VTK_VTR:
493: DMDAVTKWriteAll_VTR(dm,viewer);
494: break;
495: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]);
496: }
497: return(0);
498: }