1 #include <petsc/private/dmdaimpl.h> 2 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMDAVTKWriteAll_VTS" 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; 27 28 PetscFunctionBegin; 29 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 30 #if defined(PETSC_USE_COMPLEX) 31 SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported"); 32 #endif 33 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 34 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 35 ierr = DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);CHKERRQ(ierr); 36 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 37 ierr = DMDAGetBoundingBox(da,gmin,gmax);CHKERRQ(ierr); 38 ierr = DMGetCoordinates(da,&Coords);CHKERRQ(ierr); 39 if (Coords) { 40 PetscInt csize; 41 ierr = VecGetSize(Coords,&csize);CHKERRQ(ierr); 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 } 48 49 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 50 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 51 #if defined(PETSC_WORDS_BIGENDIAN) 52 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr); 53 #else 54 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr); 55 #endif 56 ierr = PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);CHKERRQ(ierr); 57 58 if (!rank) {ierr = PetscMalloc1(size*6,&grloc);CHKERRQ(ierr);} 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 ierr = MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);CHKERRQ(ierr); 66 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 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);CHKERRQ(ierr); 86 break; 87 case 2: 88 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm,ys+ym-1,xs,xs+xm-1,0,0);CHKERRQ(ierr); 89 break; 90 case 3: 91 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);CHKERRQ(ierr); 92 break; 93 default: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"No support for dimension %D",dim); 94 } 95 #endif 96 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);CHKERRQ(ierr); 97 ierr = PetscFPrintf(comm,fp," <Points>\n");CHKERRQ(ierr); 98 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 99 boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int); 100 ierr = PetscFPrintf(comm,fp," </Points>\n");CHKERRQ(ierr); 101 102 ierr = PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");CHKERRQ(ierr); 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 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 108 } 109 for (i=0; i<bs; i++) { 110 char buf[256]; 111 const char *fieldname; 112 ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr); 113 if (!fieldname) { 114 ierr = PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);CHKERRQ(ierr); 115 fieldname = buf; 116 } 117 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 118 boffset += nnodes*sizeof(PetscScalar) + sizeof(int); 119 } 120 } 121 ierr = PetscFPrintf(comm,fp," </PointData>\n");CHKERRQ(ierr); 122 ierr = PetscFPrintf(comm,fp," </Piece>\n");CHKERRQ(ierr); 123 } 124 ierr = PetscFPrintf(comm,fp," </StructuredGrid>\n");CHKERRQ(ierr); 125 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 126 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 127 128 /* Now write the arrays. */ 129 tag = ((PetscObject)viewer)->tag; 130 ierr = PetscMalloc2(maxnnodes*PetscMax(3,bs),&array,maxnnodes*3,&array2);CHKERRQ(ierr); 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 } 145 146 /* Write the coordinates */ 147 if (Coords) { 148 const PetscScalar *coords; 149 ierr = VecGetArrayRead(Coords,&coords);CHKERRQ(ierr); 150 if (!rank) { 151 if (r) { 152 PetscMPIInt nn; 153 ierr = MPI_Recv(array,nnodes*cdim,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 154 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 155 if (nn != nnodes*cdim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch"); 156 } else { 157 ierr = PetscMemcpy(array,coords,nnodes*cdim*sizeof(PetscScalar));CHKERRQ(ierr); 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 ierr = MPI_Send((void*)coords,nnodes*cdim,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 172 } 173 ierr = VecRestoreArrayRead(Coords,&coords);CHKERRQ(ierr); 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 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes*3,PETSC_SCALAR);CHKERRQ(ierr); 187 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; 192 193 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 194 if (!rank) { 195 if (r) { 196 PetscMPIInt nn; 197 ierr = MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 198 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 199 if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r); 200 } else { 201 ierr = PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));CHKERRQ(ierr); 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 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);CHKERRQ(ierr); 214 } 215 } else if (r == rank) { 216 ierr = MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 217 } 218 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 219 } 220 } 221 ierr = PetscFree2(array,array2);CHKERRQ(ierr); 222 ierr = PetscFree(grloc);CHKERRQ(ierr); 223 224 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 225 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 226 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "DMDAVTKWriteAll_VTR" 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; 253 254 PetscFunctionBegin; 255 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 256 #if defined(PETSC_USE_COMPLEX) 257 SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported"); 258 #endif 259 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 260 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 261 ierr = DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);CHKERRQ(ierr); 262 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 263 ierr = DMDAGetBoundingBox(da,gmin,gmax);CHKERRQ(ierr); 264 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 265 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 266 #if defined(PETSC_WORDS_BIGENDIAN) 267 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr); 268 #else 269 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr); 270 #endif 271 ierr = PetscFPrintf(comm,fp," <RectilinearGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);CHKERRQ(ierr); 272 273 if (!rank) {ierr = PetscMalloc1(size*6,&grloc);CHKERRQ(ierr);} 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 ierr = MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);CHKERRQ(ierr); 281 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 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);CHKERRQ(ierr); 298 ierr = PetscFPrintf(comm,fp," <Coordinates>\n");CHKERRQ(ierr); 299 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Xcoord\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 300 boffset += xm*sizeof(PetscScalar) + sizeof(int); 301 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Ycoord\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 302 boffset += ym*sizeof(PetscScalar) + sizeof(int); 303 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Zcoord\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 304 boffset += zm*sizeof(PetscScalar) + sizeof(int); 305 ierr = PetscFPrintf(comm,fp," </Coordinates>\n");CHKERRQ(ierr); 306 ierr = PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");CHKERRQ(ierr); 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 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 312 } 313 for (i=0; i<bs; i++) { 314 char buf[256]; 315 const char *fieldname; 316 ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr); 317 if (!fieldname) { 318 ierr = PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);CHKERRQ(ierr); 319 fieldname = buf; 320 } 321 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 322 boffset += nnodes*sizeof(PetscScalar) + sizeof(int); 323 } 324 } 325 ierr = PetscFPrintf(comm,fp," </PointData>\n");CHKERRQ(ierr); 326 ierr = PetscFPrintf(comm,fp," </Piece>\n");CHKERRQ(ierr); 327 } 328 ierr = PetscFPrintf(comm,fp," </RectilinearGrid>\n");CHKERRQ(ierr); 329 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 330 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 331 332 /* Now write the arrays. */ 333 tag = ((PetscObject)viewer)->tag; 334 ierr = PetscMalloc2(PetscMax(maxnnodes*bs,info.xm+info.ym+info.zm),&array,PetscMax(maxnnodes*3,info.xm+info.ym+info.zm),&array2);CHKERRQ(ierr); 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 ierr = DMGetCoordinates(da,&Coords);CHKERRQ(ierr); 352 if (Coords) { 353 const PetscScalar *coords; 354 ierr = VecGetArrayRead(Coords,&coords);CHKERRQ(ierr); 355 if (!rank) { 356 if (r) { 357 PetscMPIInt nn; 358 ierr = MPI_Recv(array,xm+ym+zm,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 359 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 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 ierr = MPI_Send((void*)array2,xm+ym+zm,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 395 } 396 ierr = VecRestoreArrayRead(Coords,&coords);CHKERRQ(ierr); 397 } else { /* Fabricate some coordinates using grid index */ 398 for (i=0; i<xm; i++) { 399 array[i] = xs+i; 400 } 401 for (j=0; j<ym; j++) { 402 array[j+xm] = ys+j; 403 } 404 for (k=0; k<zm; k++) { 405 array[k+xm+ym] = zs+k; 406 } 407 } 408 if (!rank) { 409 ierr = PetscViewerVTKFWrite(viewer,fp,&(array[0]),xm,PETSC_SCALAR);CHKERRQ(ierr); 410 ierr = PetscViewerVTKFWrite(viewer,fp,&(array[xm]),ym,PETSC_SCALAR);CHKERRQ(ierr); 411 ierr = PetscViewerVTKFWrite(viewer,fp,&(array[xm+ym]),zm,PETSC_SCALAR);CHKERRQ(ierr); 412 } 413 } 414 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; 419 420 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 421 if (!rank) { 422 if (r) { 423 PetscMPIInt nn; 424 ierr = MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 425 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 426 if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r); 427 } else { 428 ierr = PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));CHKERRQ(ierr); 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 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);CHKERRQ(ierr); 441 } 442 } else if (r == rank) { 443 ierr = MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 444 } 445 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 446 } 447 } 448 ierr = PetscFree2(array,array2);CHKERRQ(ierr); 449 ierr = PetscFree(grloc);CHKERRQ(ierr); 450 451 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 452 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 453 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 454 PetscFunctionReturn(0); 455 } 456 457 #undef __FUNCT__ 458 #define __FUNCT__ "DMDAVTKWriteAll" 459 /*@C 460 DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer 461 462 Collective 463 464 Input Arguments: 465 odm - DM specifying the grid layout, passed as a PetscObject 466 viewer - viewer of type VTK 467 468 Level: developer 469 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. 474 475 .seealso: PETSCVIEWERVTK 476 @*/ 477 PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer) 478 { 479 DM dm = (DM)odm; 480 PetscBool isvtk; 481 PetscErrorCode ierr; 482 483 PetscFunctionBegin; 484 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 485 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 486 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 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 ierr = DMDAVTKWriteAll_VTS(dm,viewer);CHKERRQ(ierr); 491 break; 492 case PETSC_VIEWER_VTK_VTR: 493 ierr = DMDAVTKWriteAll_VTR(dm,viewer);CHKERRQ(ierr); 494 break; 495 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]); 496 } 497 PetscFunctionReturn(0); 498 } 499