Actual source code: grvtk.c
petsc-3.4.5 2014-06-29
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: PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data;
17: PetscViewerVTKObjectLink link;
18: FILE *fp;
19: PetscMPIInt rank,size,tag;
20: DMDALocalInfo info;
21: PetscInt dim,mx,my,mz,bs,boffset,maxnnodes,i,j,k,f,r;
22: PetscInt rloc[6],(*grloc)[6] = NULL;
23: PetscScalar *array,*array2;
24: PetscReal gmin[3],gmax[3];
25: PetscErrorCode ierr;
28: PetscObjectGetComm((PetscObject)da,&comm);
29: #if defined(PETSC_USE_COMPLEX)
30: SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported");
31: #endif
32: MPI_Comm_size(comm,&size);
33: MPI_Comm_rank(comm,&rank);
34: DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);
35: DMDAGetLocalInfo(da,&info);
36: DMDAGetBoundingBox(da,gmin,gmax);
38: PetscFOpen(comm,vtk->filename,"wb",&fp);
39: PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
40: #if defined(PETSC_WORDS_BIGENDIAN)
41: PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
42: #else
43: PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
44: #endif
45: PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);
47: if (!rank) {PetscMalloc(size*6*sizeof(PetscInt),&grloc);}
48: rloc[0] = info.xs;
49: rloc[1] = info.xm;
50: rloc[2] = info.ys;
51: rloc[3] = info.ym;
52: rloc[4] = info.zs;
53: rloc[5] = info.zm;
54: MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);
56: /* Write XML header */
57: maxnnodes = 0; /* Used for the temporary array size on rank 0 */
58: boffset = 0; /* Offset into binary file */
59: for (r=0; r<size; r++) {
60: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
61: if (!rank) {
62: xs = grloc[r][0];
63: xm = grloc[r][1];
64: ys = grloc[r][2];
65: ym = grloc[r][3];
66: zs = grloc[r][4];
67: zm = grloc[r][5];
68: nnodes = xm*ym*zm;
69: }
70: maxnnodes = PetscMax(maxnnodes,nnodes);
71: #if 0
72: switch (dim) {
73: case 1:
74: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);
75: break;
76: case 2:
77: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm,ys+ym-1,xs,xs+xm-1,0,0);
78: break;
79: case 3:
80: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
81: break;
82: default: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"No support for dimension %D",dim);
83: }
84: #endif
85: PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
86: PetscFPrintf(comm,fp," <Points>\n");
87: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
88: boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int);
89: PetscFPrintf(comm,fp," </Points>\n");
91: PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");
92: for (link=vtk->link; link; link=link->next) {
93: Vec X = (Vec)link->vec;
94: const char *vecname = "";
95: 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. */
96: PetscObjectGetName((PetscObject)X,&vecname);
97: }
98: for (i=0; i<bs; i++) {
99: char buf[256];
100: const char *fieldname;
101: DMDAGetFieldName(da,i,&fieldname);
102: if (!fieldname) {
103: PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);
104: fieldname = buf;
105: }
106: PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
107: boffset += nnodes*sizeof(PetscScalar) + sizeof(int);
108: }
109: }
110: PetscFPrintf(comm,fp," </PointData>\n");
111: PetscFPrintf(comm,fp," </Piece>\n");
112: }
113: PetscFPrintf(comm,fp," </StructuredGrid>\n");
114: PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");
115: PetscFPrintf(comm,fp,"_");
117: /* Now write the arrays. */
118: tag = ((PetscObject)viewer)->tag;
119: PetscMalloc2(maxnnodes*PetscMax(3,bs),PetscScalar,&array,maxnnodes*3,PetscScalar,&array2);
120: for (r=0; r<size; r++) {
121: MPI_Status status;
122: PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
123: if (!rank) {
124: xs = grloc[r][0];
125: xm = grloc[r][1];
126: ys = grloc[r][2];
127: ym = grloc[r][3];
128: zs = grloc[r][4];
129: zm = grloc[r][5];
130: nnodes = xm*ym*zm;
131: } else if (r == rank) {
132: nnodes = info.xm*info.ym*info.zm;
133: }
135: { /* Write the coordinates */
136: Vec Coords;
137: DMGetCoordinates(da,&Coords);
138: if (Coords) {
139: const PetscScalar *coords;
140: VecGetArrayRead(Coords,&coords);
141: if (!rank) {
142: if (r) {
143: PetscMPIInt nn;
144: MPI_Recv(array,nnodes*dim,MPIU_SCALAR,r,tag,comm,&status);
145: MPI_Get_count(&status,MPIU_SCALAR,&nn);
146: if (nn != nnodes*dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
147: } else {
148: PetscMemcpy(array,coords,nnodes*dim*sizeof(PetscScalar));
149: }
150: /* Transpose coordinates to VTK (C-style) ordering */
151: for (k=0; k<zm; k++) {
152: for (j=0; j<ym; j++) {
153: for (i=0; i<xm; i++) {
154: PetscInt Iloc = i+xm*(j+ym*k);
155: array2[Iloc*3+0] = array[Iloc*dim + 0];
156: array2[Iloc*3+1] = dim > 1 ? array[Iloc*dim + 1] : 0;
157: array2[Iloc*3+2] = dim > 2 ? array[Iloc*dim + 2] : 0;
158: }
159: }
160: }
161: } else if (r == rank) {
162: MPI_Send((void*)coords,nnodes*dim,MPIU_SCALAR,0,tag,comm);
163: }
164: VecRestoreArrayRead(Coords,&coords);
165: } else { /* Fabricate some coordinates using grid index */
166: for (k=0; k<zm; k++) {
167: for (j=0; j<ym; j++) {
168: for (i=0; i<xm; i++) {
169: PetscInt Iloc = i+xm*(j+ym*k);
170: array2[Iloc*3+0] = xs+i;
171: array2[Iloc*3+1] = ys+j;
172: array2[Iloc*3+2] = zs+k;
173: }
174: }
175: }
176: }
177: PetscViewerVTKFWrite(viewer,fp,array2,nnodes*3,PETSC_SCALAR);
178: }
180: /* Write each of the objects queued up for this file */
181: for (link=vtk->link; link; link=link->next) {
182: Vec X = (Vec)link->vec;
183: const PetscScalar *x;
185: VecGetArrayRead(X,&x);
186: if (!rank) {
187: if (r) {
188: PetscMPIInt nn;
189: MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);
190: MPI_Get_count(&status,MPIU_SCALAR,&nn);
191: if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r);
192: } else {
193: PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));
194: }
195: for (f=0; f<bs; f++) {
196: /* Extract and transpose the f'th field */
197: for (k=0; k<zm; k++) {
198: for (j=0; j<ym; j++) {
199: for (i=0; i<xm; i++) {
200: PetscInt Iloc = i+xm*(j+ym*k);
201: array2[Iloc] = array[Iloc*bs + f];
202: }
203: }
204: }
205: PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);
206: }
207: } else if (r == rank) {
208: MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);
209: }
210: VecRestoreArrayRead(X,&x);
211: }
212: }
213: PetscFree2(array,array2);
214: PetscFree(grloc);
216: PetscFPrintf(comm,fp,"\n </AppendedData>\n");
217: PetscFPrintf(comm,fp,"</VTKFile>\n");
218: PetscFClose(comm,fp);
219: return(0);
220: }
225: /*@C
226: DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
228: Collective
230: Input Arguments:
231: odm - DM specifying the grid layout, passed as a PetscObject
232: viewer - viewer of type VTK
234: Level: developer
236: Note:
237: This function is a callback used by the VTK viewer to actually write the file.
238: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
239: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
241: .seealso: PETSCVIEWERVTK
242: @*/
243: PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer)
244: {
245: DM dm = (DM)odm;
246: PetscBool isvtk;
252: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);
253: if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name);
254: switch (viewer->format) {
255: case PETSC_VIEWER_VTK_VTS:
256: DMDAVTKWriteAll_VTS(dm,viewer);
257: break;
258: default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]);
259: }
260: return(0);
261: }