1 #include <petsc/private/dmdaimpl.h> 2 /* 3 Note that the API for using PETSCVIEWERVTK is totally wrong since its use requires 4 including the private vtkvimpl.h file. The code should be refactored. 5 */ 6 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h> 7 8 static PetscErrorCode DMDAVTKWriteAll_VTS(DM da,PetscViewer viewer) 9 { 10 #if defined(PETSC_USE_REAL_SINGLE) 11 const char precision[] = "Float32"; 12 #elif defined(PETSC_USE_REAL_DOUBLE) 13 const char precision[] = "Float64"; 14 #else 15 const char precision[] = "UnknownPrecision"; 16 #endif 17 MPI_Comm comm; 18 Vec Coords; 19 PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data; 20 PetscViewerVTKObjectLink link; 21 FILE *fp; 22 PetscMPIInt rank,size,tag; 23 DMDALocalInfo info; 24 PetscInt dim,mx,my,mz,cdim,bs,boffset,maxnnodes,i,j,k,f,r; 25 PetscInt rloc[6],(*grloc)[6] = NULL; 26 PetscScalar *array,*array2; 27 PetscReal gmin[3],gmax[3]; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 32 #if defined(PETSC_USE_COMPLEX) 33 SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported"); 34 #endif 35 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 36 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 37 ierr = DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);CHKERRQ(ierr); 38 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 39 ierr = DMDAGetBoundingBox(da,gmin,gmax);CHKERRQ(ierr); 40 ierr = DMGetCoordinates(da,&Coords);CHKERRQ(ierr); 41 if (Coords) { 42 PetscInt csize; 43 ierr = VecGetSize(Coords,&csize);CHKERRQ(ierr); 44 if (csize % (mx*my*mz)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate vector size mismatch"); 45 cdim = csize/(mx*my*mz); 46 if (cdim < dim || cdim > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate vector size mismatch"); 47 } else { 48 cdim = dim; 49 } 50 51 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 52 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 53 #if defined(PETSC_WORDS_BIGENDIAN) 54 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr); 55 #else 56 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr); 57 #endif 58 ierr = PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);CHKERRQ(ierr); 59 60 if (!rank) {ierr = PetscMalloc1(size*6,&grloc);CHKERRQ(ierr);} 61 rloc[0] = info.xs; 62 rloc[1] = info.xm; 63 rloc[2] = info.ys; 64 rloc[3] = info.ym; 65 rloc[4] = info.zs; 66 rloc[5] = info.zm; 67 ierr = MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);CHKERRQ(ierr); 68 69 /* Write XML header */ 70 maxnnodes = 0; /* Used for the temporary array size on rank 0 */ 71 boffset = 0; /* Offset into binary file */ 72 for (r=0; r<size; r++) { 73 PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0; 74 if (!rank) { 75 xs = grloc[r][0]; 76 xm = grloc[r][1]; 77 ys = grloc[r][2]; 78 ym = grloc[r][3]; 79 zs = grloc[r][4]; 80 zm = grloc[r][5]; 81 nnodes = xm*ym*zm; 82 } 83 maxnnodes = PetscMax(maxnnodes,nnodes); 84 #if 0 85 switch (dim) { 86 case 1: 87 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);CHKERRQ(ierr); 88 break; 89 case 2: 90 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); 91 break; 92 case 3: 93 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); 94 break; 95 default: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"No support for dimension %D",dim); 96 } 97 #endif 98 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); 99 ierr = PetscFPrintf(comm,fp," <Points>\n");CHKERRQ(ierr); 100 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 101 boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int); 102 ierr = PetscFPrintf(comm,fp," </Points>\n");CHKERRQ(ierr); 103 104 ierr = PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");CHKERRQ(ierr); 105 for (link=vtk->link; link; link=link->next) { 106 Vec X = (Vec)link->vec; 107 const char *vecname = ""; 108 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. */ 109 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 110 } 111 for (i=0; i<bs; i++) { 112 char buf[256]; 113 const char *fieldname; 114 ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr); 115 if (!fieldname) { 116 ierr = PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);CHKERRQ(ierr); 117 fieldname = buf; 118 } 119 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 120 boffset += nnodes*sizeof(PetscScalar) + sizeof(int); 121 } 122 } 123 ierr = PetscFPrintf(comm,fp," </PointData>\n");CHKERRQ(ierr); 124 ierr = PetscFPrintf(comm,fp," </Piece>\n");CHKERRQ(ierr); 125 } 126 ierr = PetscFPrintf(comm,fp," </StructuredGrid>\n");CHKERRQ(ierr); 127 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 128 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 129 130 /* Now write the arrays. */ 131 tag = ((PetscObject)viewer)->tag; 132 ierr = PetscMalloc2(maxnnodes*PetscMax(3,bs),&array,maxnnodes*3,&array2);CHKERRQ(ierr); 133 for (r=0; r<size; r++) { 134 MPI_Status status; 135 PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0; 136 if (!rank) { 137 xs = grloc[r][0]; 138 xm = grloc[r][1]; 139 ys = grloc[r][2]; 140 ym = grloc[r][3]; 141 zs = grloc[r][4]; 142 zm = grloc[r][5]; 143 nnodes = xm*ym*zm; 144 } else if (r == rank) { 145 nnodes = info.xm*info.ym*info.zm; 146 } 147 148 /* Write the coordinates */ 149 if (Coords) { 150 const PetscScalar *coords; 151 ierr = VecGetArrayRead(Coords,&coords);CHKERRQ(ierr); 152 if (!rank) { 153 if (r) { 154 PetscMPIInt nn; 155 ierr = MPI_Recv(array,nnodes*cdim,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 156 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 157 if (nn != nnodes*cdim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch"); 158 } else { 159 ierr = PetscMemcpy(array,coords,nnodes*cdim*sizeof(PetscScalar));CHKERRQ(ierr); 160 } 161 /* Transpose coordinates to VTK (C-style) ordering */ 162 for (k=0; k<zm; k++) { 163 for (j=0; j<ym; j++) { 164 for (i=0; i<xm; i++) { 165 PetscInt Iloc = i+xm*(j+ym*k); 166 array2[Iloc*3+0] = array[Iloc*cdim + 0]; 167 array2[Iloc*3+1] = cdim > 1 ? array[Iloc*cdim + 1] : 0.0; 168 array2[Iloc*3+2] = cdim > 2 ? array[Iloc*cdim + 2] : 0.0; 169 } 170 } 171 } 172 } else if (r == rank) { 173 ierr = MPI_Send((void*)coords,nnodes*cdim,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 174 } 175 ierr = VecRestoreArrayRead(Coords,&coords);CHKERRQ(ierr); 176 } else { /* Fabricate some coordinates using grid index */ 177 for (k=0; k<zm; k++) { 178 for (j=0; j<ym; j++) { 179 for (i=0; i<xm; i++) { 180 PetscInt Iloc = i+xm*(j+ym*k); 181 array2[Iloc*3+0] = xs+i; 182 array2[Iloc*3+1] = ys+j; 183 array2[Iloc*3+2] = zs+k; 184 } 185 } 186 } 187 } 188 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes*3,MPIU_SCALAR);CHKERRQ(ierr); 189 190 /* Write each of the objects queued up for this file */ 191 for (link=vtk->link; link; link=link->next) { 192 Vec X = (Vec)link->vec; 193 const PetscScalar *x; 194 195 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 196 if (!rank) { 197 if (r) { 198 PetscMPIInt nn; 199 ierr = MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 200 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 201 if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r); 202 } else { 203 ierr = PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));CHKERRQ(ierr); 204 } 205 for (f=0; f<bs; f++) { 206 /* Extract and transpose the f'th field */ 207 for (k=0; k<zm; k++) { 208 for (j=0; j<ym; j++) { 209 for (i=0; i<xm; i++) { 210 PetscInt Iloc = i+xm*(j+ym*k); 211 array2[Iloc] = array[Iloc*bs + f]; 212 } 213 } 214 } 215 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes,MPIU_SCALAR);CHKERRQ(ierr); 216 } 217 } else if (r == rank) { 218 ierr = MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 219 } 220 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 221 } 222 } 223 ierr = PetscFree2(array,array2);CHKERRQ(ierr); 224 ierr = PetscFree(grloc);CHKERRQ(ierr); 225 226 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 227 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 228 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 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,MPIU_SCALAR);CHKERRQ(ierr); 410 ierr = PetscViewerVTKFWrite(viewer,fp,&(array[xm]),ym,MPIU_SCALAR);CHKERRQ(ierr); 411 ierr = PetscViewerVTKFWrite(viewer,fp,&(array[xm+ym]),zm,MPIU_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,MPIU_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 /*@C 458 DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer 459 460 Collective 461 462 Input Arguments: 463 odm - DM specifying the grid layout, passed as a PetscObject 464 viewer - viewer of type VTK 465 466 Level: developer 467 468 Note: 469 This function is a callback used by the VTK viewer to actually write the file. 470 The reason for this odd model is that the VTK file format does not provide any way to write one field at a time. 471 Instead, metadata for the entire file needs to be available up-front before you can start writing the file. 472 473 .seealso: PETSCVIEWERVTK 474 @*/ 475 PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer) 476 { 477 DM dm = (DM)odm; 478 PetscBool isvtk; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 483 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 484 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 485 if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name); 486 switch (viewer->format) { 487 case PETSC_VIEWER_VTK_VTS: 488 ierr = DMDAVTKWriteAll_VTS(dm,viewer);CHKERRQ(ierr); 489 break; 490 case PETSC_VIEWER_VTK_VTR: 491 ierr = DMDAVTKWriteAll_VTR(dm,viewer);CHKERRQ(ierr); 492 break; 493 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]); 494 } 495 PetscFunctionReturn(0); 496 } 497