1 #include <petsc-private/daimpl.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 = ((PetscObject)da)->comm; 16 PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data; 17 PetscViewerVTKObjectLink link; 18 FILE *fp; 19 PetscMPIInt rank,size,tag; 20 DMDALocalInfo info; 21 PetscInt dim,mx,my,mz,bs,boffset,maxnnodes,i,j,k,f,r; 22 PetscInt rloc[6],(*grloc)[6] = PETSC_NULL; 23 PetscScalar *array,*array2; 24 PetscReal gmin[3],gmax[3]; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 #if defined(PETSC_USE_COMPLEX) 29 SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported"); 30 #endif 31 32 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 33 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 34 ierr = DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);CHKERRQ(ierr); 35 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 36 ierr = DMDAGetBoundingBox(da,gmin,gmax);CHKERRQ(ierr); 37 38 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 39 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 40 #if defined(PETSC_WORDS_BIGENDIAN) 41 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr); 42 #else 43 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr); 44 #endif 45 ierr = PetscFPrintf(comm,fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);CHKERRQ(ierr); 46 47 if (!rank) {ierr = PetscMalloc(size*6*sizeof(PetscInt),&grloc);CHKERRQ(ierr);} 48 rloc[0] = info.xs; 49 rloc[1] = info.xm; 50 rloc[2] = info.ys; 51 rloc[3] = info.ym; 52 rloc[4] = info.zs; 53 rloc[5] = info.zm; 54 ierr = MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);CHKERRQ(ierr); 55 56 /* Write XML header */ 57 maxnnodes = 0; /* Used for the temporary array size on rank 0 */ 58 boffset = 0; /* Offset into binary file */ 59 for (r=0; r<size; r++) { 60 PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0; 61 if (!rank) { 62 xs = grloc[r][0]; 63 xm = grloc[r][1]; 64 ys = grloc[r][2]; 65 ym = grloc[r][3]; 66 zs = grloc[r][4]; 67 zm = grloc[r][5]; 68 nnodes = xm*ym*zm; 69 } 70 maxnnodes = PetscMax(maxnnodes,nnodes); 71 #if 0 72 switch (dim) { 73 case 1: 74 ierr = PetscFPrintf(comm,fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);CHKERRQ(ierr); 75 break; 76 case 2: 77 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); 78 break; 79 case 3: 80 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); 81 break; 82 default: SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"No support for dimension %D",dim); 83 } 84 #endif 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 const char *vecname = ""; 95 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. */ 96 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 97 } 98 for (i=0; i<bs; i++) { 99 char buf[256]; 100 const char *fieldname; 101 ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr); 102 if (!fieldname) { 103 ierr = PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);CHKERRQ(ierr); 104 fieldname = buf; 105 } 106 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 107 boffset += nnodes*sizeof(PetscScalar) + sizeof(int); 108 } 109 } 110 ierr = PetscFPrintf(comm,fp," </PointData>\n");CHKERRQ(ierr); 111 ierr = PetscFPrintf(comm,fp," </Piece>\n");CHKERRQ(ierr); 112 } 113 ierr = PetscFPrintf(comm,fp," </StructuredGrid>\n");CHKERRQ(ierr); 114 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 115 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 116 117 /* Now write the arrays. */ 118 tag = ((PetscObject)viewer)->tag; 119 ierr = PetscMalloc2(maxnnodes*PetscMax(3,bs),PetscScalar,&array,maxnnodes*3,PetscScalar,&array2);CHKERRQ(ierr); 120 for (r=0; r<size; r++) { 121 MPI_Status status; 122 PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0; 123 if (!rank) { 124 xs = grloc[r][0]; 125 xm = grloc[r][1]; 126 ys = grloc[r][2]; 127 ym = grloc[r][3]; 128 zs = grloc[r][4]; 129 zm = grloc[r][5]; 130 nnodes = xm*ym*zm; 131 } else if (r == rank) { 132 nnodes = info.xm*info.ym*info.zm; 133 } 134 135 { /* Write the coordinates */ 136 Vec Coords; 137 ierr = DMGetCoordinates(da,&Coords);CHKERRQ(ierr); 138 if (Coords) { 139 const PetscScalar *coords; 140 ierr = VecGetArrayRead(Coords,&coords);CHKERRQ(ierr); 141 if (!rank) { 142 if (r) { 143 PetscMPIInt nn; 144 ierr = MPI_Recv(array,nnodes*dim,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 145 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 146 if (nn != nnodes*dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch"); 147 } else { 148 ierr = PetscMemcpy(array,coords,nnodes*dim*sizeof(PetscScalar));CHKERRQ(ierr); 149 } 150 /* Transpose coordinates to VTK (C-style) ordering */ 151 for (k=0; k<zm; k++) { 152 for (j=0; j<ym; j++) { 153 for (i=0; i<xm; i++) { 154 PetscInt Iloc = i+xm*(j+ym*k); 155 array2[Iloc*3+0] = array[Iloc*dim + 0]; 156 array2[Iloc*3+1] = dim > 1 ? array[Iloc*dim + 1] : 0; 157 array2[Iloc*3+2] = dim > 2 ? array[Iloc*dim + 2] : 0; 158 } 159 } 160 } 161 } else if (r == rank) { 162 ierr = MPI_Send((void*)coords,nnodes*dim,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 163 } 164 ierr = VecRestoreArrayRead(Coords,&coords);CHKERRQ(ierr); 165 } else { /* Fabricate some coordinates using grid index */ 166 for (k=0; k<zm; k++) { 167 for (j=0; j<ym; j++) { 168 for (i=0; i<xm; i++) { 169 PetscInt Iloc = i+xm*(j+ym*k); 170 array2[Iloc*3+0] = xs+i; 171 array2[Iloc*3+1] = ys+j; 172 array2[Iloc*3+2] = zs+k; 173 } 174 } 175 } 176 } 177 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes*3,PETSC_SCALAR);CHKERRQ(ierr); 178 } 179 180 /* Write each of the objects queued up for this file */ 181 for (link=vtk->link; link; link=link->next) { 182 Vec X = (Vec)link->vec; 183 const PetscScalar *x; 184 185 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 186 if (!rank) { 187 if (r) { 188 PetscMPIInt nn; 189 ierr = MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);CHKERRQ(ierr); 190 ierr = MPI_Get_count(&status,MPIU_SCALAR,&nn);CHKERRQ(ierr); 191 if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r); 192 } else { 193 ierr = PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));CHKERRQ(ierr); 194 } 195 for (f=0; f<bs; f++) { 196 /* Extract and transpose the f'th field */ 197 for (k=0; k<zm; k++) { 198 for (j=0; j<ym; j++) { 199 for (i=0; i<xm; i++) { 200 PetscInt Iloc = i+xm*(j+ym*k); 201 array2[Iloc] = array[Iloc*bs + f]; 202 } 203 } 204 } 205 ierr = PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);CHKERRQ(ierr); 206 } 207 } else if (r == rank) { 208 ierr = MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 209 } 210 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 211 } 212 } 213 ierr = PetscFree2(array,array2);CHKERRQ(ierr); 214 ierr = PetscFree(grloc);CHKERRQ(ierr); 215 216 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 217 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 218 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "DMDAVTKWriteAll" 225 /*@C 226 DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer 227 228 Collective 229 230 Input Arguments: 231 odm - DM specifying the grid layout, passed as a PetscObject 232 viewer - viewer of type VTK 233 234 Level: developer 235 236 Note: 237 This function is a callback used by the VTK viewer to actually write the file. 238 The reason for this odd model is that the VTK file format does not provide any way to write one field at a time. 239 Instead, metadata for the entire file needs to be available up-front before you can start writing the file. 240 241 .seealso: PETSCVIEWERVTK 242 @*/ 243 PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer) 244 { 245 DM dm = (DM)odm; 246 PetscBool isvtk; 247 PetscErrorCode ierr; 248 249 PetscFunctionBegin; 250 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 251 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 252 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr); 253 if (!isvtk) SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name); 254 switch (viewer->format) { 255 case PETSC_VIEWER_VTK_VTS: 256 ierr = DMDAVTKWriteAll_VTS(dm,viewer);CHKERRQ(ierr); 257 break; 258 default: SETERRQ1(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]); 259 } 260 PetscFunctionReturn(0); 261 } 262