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