1 #include <petsc-private/daimpl.h> 2 #include <../src/sys/viewer/impls/vtk/vtkvimpl.h> 3 4 #if defined(PETSC_HAVE_STDINT_H) /* The VTK format requires a 32-bit integer */ 5 typedef int32_t PetscVTKInt; 6 #else 7 typedef int PetscVTKInt; 8 #endif 9 #define PETSC_VTK_INT_MAX 2147483647 10 #define PETSC_VTK_INT_MIN -2147483647 11 #if defined(PETSC_USE_64BIT_INDICES) 12 # define PetscVTKIntCheck(a) if ((a) > PETSC_VTK_INT_MAX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Array too long for 32-bit VTK binary format") 13 # define PetscVTKIntCast(a) (PetscVTKInt)(a);PetscVTKIntCheck(a) 14 #else 15 # define PetscVTKIntCheck(a) 16 # define PetscVTKIntCast(a) a 17 #endif 18 19 #undef __FUNCT__ 20 #define __FUNCT__ "PetscFWrite_VTK" 21 /* Write binary data preceded by 32-bit int length (in bytes), does not do byte swapping. */ 22 static PetscErrorCode PetscFWrite_VTK(MPI_Comm comm,FILE *fp,void *data,PetscInt n,PetscDataType dtype) 23 { 24 PetscErrorCode ierr; 25 PetscMPIInt rank; 26 27 PetscFunctionBegin; 28 if (n < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Trying to write a negative amount of data %D",n); 29 if (!n) PetscFunctionReturn(0); 30 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 31 if (!rank) { 32 size_t count; 33 PetscInt size; 34 PetscVTKInt bytes; 35 switch (dtype) { 36 case PETSC_DOUBLE: 37 size = sizeof(double); 38 break; 39 case PETSC_FLOAT: 40 size = sizeof(float); 41 break; 42 case PETSC_INT: 43 size = sizeof(PetscInt); 44 break; 45 case PETSC_CHAR: 46 size = sizeof(char); 47 break; 48 default: SETERRQ(comm,PETSC_ERR_SUP,"Data type not supported"); 49 } 50 bytes = PetscVTKIntCast(size*n); 51 52 count = fwrite(&bytes,sizeof(int),1,fp); 53 if (count != 1) { 54 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_WRITE,"Error writing byte count"); 55 } 56 count = fwrite(data,size,(size_t)n,fp); 57 if ((PetscInt)count != n) { 58 SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_FILE_WRITE,"Wrote %D/%D array members of size %D",(PetscInt)count,n,(PetscInt)size); 59 } 60 } 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "DMDAVTKWriteAll_VTS" 66 static PetscErrorCode DMDAVTKWriteAll_VTS(DM da,PetscViewer viewer) 67 { 68 #if defined(PETSC_USE_REAL_SINGLE) 69 const char precision[] = "Float32"; 70 #elif defined(PETSC_USE_REAL_DOUBLE) 71 const char precision[] = "Float64"; 72 #else 73 const char precision[] = "UnknownPrecision"; 74 #endif 75 MPI_Comm comm = ((PetscObject)da)->comm; 76 PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data; 77 PetscViewerVTKObjectLink link; 78 FILE *fp; 79 PetscMPIInt rank,size,tag; 80 DMDALocalInfo info; 81 PetscInt dim,mx,my,mz,bs,boffset,maxnnodes,i,j,k,f,r; 82 PetscInt rloc[6],(*grloc)[6] = PETSC_NULL; 83 PetscScalar *array,*array2; 84 PetscReal gmin[3],gmax[3]; 85 PetscErrorCode ierr; 86 87 PetscFunctionBegin; 88 #if defined(PETSC_USE_COMPLEX) 89 SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP,"Complex values not supported"); 90 #endif 91 92 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 93 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 94 ierr = DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);CHKERRQ(ierr); 95 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 96 ierr = DMDAGetBoundingBox(da,gmin,gmax);CHKERRQ(ierr); 97 98 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 99 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 100 #ifdef PETSC_WORDS_BIGENDIAN 101 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr); 102 #else 103 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr); 104 #endif 105 ierr = PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);CHKERRQ(ierr); 106 107 if (!rank) {ierr = PetscMalloc(size*6*sizeof(PetscInt),&grloc);CHKERRQ(ierr);} 108 rloc[0] = info.xs; 109 rloc[1] = info.xm; 110 rloc[2] = info.ys; 111 rloc[3] = info.ym; 112 rloc[4] = info.zs; 113 rloc[5] = info.zm; 114 ierr = MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);CHKERRQ(ierr); 115 116 /* Write XML header */ 117 maxnnodes = 0; /* Used for the temporary array size on rank 0 */ 118 boffset = 0; /* Offset into binary file */ 119 for (r=0; r<size; r++) { 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 } 130 maxnnodes = PetscMax(maxnnodes,nnodes); 131 #if 0 132 switch (dim) { 133 case 1: 134 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);CHKERRQ(ierr); 135 break; 136 case 2: 137 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); 138 break; 139 case 3: 140 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); 141 break; 142 default: SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"No support for dimension %D",dim); 143 } 144 #endif 145 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); 146 ierr = PetscFPrintf(comm,fp," <Points>\n");CHKERRQ(ierr); 147 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 148 boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int); 149 ierr = PetscFPrintf(comm,fp," </Points>\n");CHKERRQ(ierr); 150 151 ierr = PetscFPrintf(comm,fp," <PointData Scalars=\"ScalarPointData\">\n");CHKERRQ(ierr); 152 for (link=vtk->link; link; link=link->next) { 153 Vec X = (Vec)link->vec; 154 const char *vecname = ""; 155 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. */ 156 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 157 } 158 for (i=0; i<bs; i++) { 159 char buf[256]; 160 const char *fieldname; 161 ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr); 162 if (!fieldname) { 163 ierr = PetscSNPrintf(buf,sizeof buf,"Unnamed%D",i);CHKERRQ(ierr); 164 fieldname = buf; 165 } 166 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 167 boffset += nnodes*sizeof(PetscScalar) + sizeof(int); 168 } 169 } 170 ierr = PetscFPrintf(comm,fp," </PointData>\n");CHKERRQ(ierr); 171 ierr = PetscFPrintf(comm,fp," </Piece>\n");CHKERRQ(ierr); 172 } 173 ierr = PetscFPrintf(comm,fp," </StructuredGrid>\n");CHKERRQ(ierr); 174 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 175 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 176 177 /* Now write the arrays. */ 178 tag = ((PetscObject)viewer)->tag; 179 ierr = PetscMalloc2(maxnnodes*PetscMax(3,bs),PetscScalar,&array,maxnnodes*3,PetscScalar,&array2);CHKERRQ(ierr); 180 for (r=0; r<size; r++) { 181 MPI_Status status; 182 PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0; 183 if (!rank) { 184 xs = grloc[r][0]; 185 xm = grloc[r][1]; 186 ys = grloc[r][2]; 187 ym = grloc[r][3]; 188 zs = grloc[r][4]; 189 zm = grloc[r][5]; 190 nnodes = xm*ym*zm; 191 } else if (r == rank) { 192 nnodes = info.xm*info.ym*info.zm; 193 } 194 195 { /* Write the coordinates */ 196 Vec Coords; 197 ierr = DMDAGetCoordinates(da,&Coords);CHKERRQ(ierr); 198 if (Coords) { 199 const PetscScalar *coords; 200 ierr = VecGetArrayRead(Coords,&coords);CHKERRQ(ierr); 201 if (!rank) { 202 if (r) { 203 PetscMPIInt nn; 204 ierr = MPI_Recv(array,nnodes*dim,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 205 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 206 if (nn != nnodes*dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch"); 207 } else { 208 ierr = PetscMemcpy(array,coords,nnodes*dim*sizeof(PetscScalar));CHKERRQ(ierr); 209 } 210 /* Transpose coordinates to VTK (C-style) ordering */ 211 for (k=0; k<zm; k++) { 212 for (j=0; j<ym; j++) { 213 for (i=0; i<xm; i++) { 214 PetscInt Iloc = i+xm*(j+ym*k); 215 array2[Iloc*3+0] = array[Iloc*dim + 0]; 216 array2[Iloc*3+1] = dim > 1 ? array[Iloc*dim + 1] : 0; 217 array2[Iloc*3+2] = dim > 2 ? array[Iloc*dim + 2] : 0; 218 } 219 } 220 } 221 } else if (r == rank) { 222 ierr = MPI_Send((void*)coords,nnodes*dim,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 223 } 224 ierr = VecRestoreArrayRead(Coords,&coords);CHKERRQ(ierr); 225 } else { /* Fabricate some coordinates using grid index */ 226 for (k=0; k<zm; k++) { 227 for (j=0; j<ym; j++) { 228 for (i=0; i<xm; i++) { 229 PetscInt Iloc = i+xm*(j+ym*k); 230 array2[Iloc*3+0] = xs+i; 231 array2[Iloc*3+1] = ys+j; 232 array2[Iloc*3+2] = zs+k; 233 } 234 } 235 } 236 } 237 ierr = PetscFWrite_VTK(comm,fp,array2,nnodes*3,PETSC_SCALAR);CHKERRQ(ierr); 238 } 239 240 /* Write each of the objects queued up for this file */ 241 for (link=vtk->link; link; link=link->next) { 242 Vec X = (Vec)link->vec; 243 const PetscScalar *x; 244 245 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 246 if (!rank) { 247 if (r) { 248 PetscMPIInt nn; 249 ierr = MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 250 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 251 if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r); 252 } else { 253 ierr = PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));CHKERRQ(ierr); 254 } 255 for (f=0; f<bs; f++) { 256 /* Extract and transpose the f'th field */ 257 for (k=0; k<zm; k++) { 258 for (j=0; j<ym; j++) { 259 for (i=0; i<xm; i++) { 260 PetscInt Iloc = i+xm*(j+ym*k); 261 array2[Iloc] = array[Iloc*bs + f]; 262 } 263 } 264 } 265 ierr = PetscFWrite_VTK(comm,fp,array2,nnodes,PETSC_SCALAR);CHKERRQ(ierr); 266 } 267 } else if (r == rank) { 268 ierr = MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 269 } 270 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 271 } 272 } 273 ierr = PetscFree2(array,array2);CHKERRQ(ierr); 274 ierr = PetscFree(grloc);CHKERRQ(ierr); 275 276 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 277 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 278 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 282 283 #undef __FUNCT__ 284 #define __FUNCT__ "DMDAVTKWriteAll" 285 /*@C 286 DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer 287 288 Collective 289 290 Input Arguments: 291 odm - DM specifying the grid layout, passed as a PetscObject 292 viewer - viewer of type VTK 293 294 Level: developer 295 296 Note: 297 This function is a callback used by the VTK viewer to actually write the file. 298 The reason for this odd model is that the VTK file format does not provide any way to write one field at a time. 299 Instead, metadata for the entire file needs to be available up-front before you can start writing the file. 300 301 .seealso: PETSCVIEWERVTK 302 @*/ 303 PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer) 304 { 305 DM dm = (DM)odm; 306 PetscBool isvtk; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 311 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 312 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 313 if (!isvtk) SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name); 314 switch (viewer->format) { 315 case PETSC_VIEWER_VTK_VTS: 316 ierr = DMDAVTKWriteAll_VTS(dm,viewer);CHKERRQ(ierr); 317 break; 318 default: SETERRQ1(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]); 319 } 320 PetscFunctionReturn(0); 321 } 322