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