1 #include <petsc/private/dmpleximpl.h> 2 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h> 3 4 typedef struct { 5 PetscInt nvertices; 6 PetscInt ncells; 7 PetscInt nconn; /* number of entries in cell->vertex connectivity array */ 8 } PieceInfo; 9 10 #if defined(PETSC_USE_REAL_SINGLE) 11 static const char precision[] = "Float32"; 12 #elif defined(PETSC_USE_REAL_DOUBLE) 13 static const char precision[] = "Float64"; 14 #else 15 static const char precision[] = "UnknownPrecision"; 16 #endif 17 18 static PetscErrorCode TransferWrite(PetscViewer viewer,FILE *fp,PetscMPIInt srank,PetscMPIInt root,const void *send,void *recv,PetscMPIInt count,PetscDataType datatype,PetscMPIInt tag) 19 { 20 PetscMPIInt rank; 21 PetscErrorCode ierr; 22 MPI_Comm comm; 23 MPI_Datatype mpidatatype; 24 25 PetscFunctionBegin; 26 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 27 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 28 ierr = PetscDataTypeToMPIDataType(datatype,&mpidatatype);CHKERRQ(ierr); 29 30 if (rank == srank && rank != root) { 31 ierr = MPI_Send((void*)send,count,mpidatatype,root,tag,comm);CHKERRQ(ierr); 32 } else if (rank == root) { 33 const void *buffer; 34 if (root == srank) { /* self */ 35 buffer = send; 36 } else { 37 MPI_Status status; 38 PetscMPIInt nrecv; 39 ierr = MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);CHKERRQ(ierr); 40 ierr = MPI_Get_count(&status,mpidatatype,&nrecv);CHKERRQ(ierr); 41 if (count != nrecv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch"); 42 buffer = recv; 43 } 44 ierr = PetscViewerVTKFWrite(viewer,fp,buffer,count,datatype);CHKERRQ(ierr); 45 } 46 PetscFunctionReturn(0); 47 } 48 49 static PetscErrorCode DMPlexGetVTKConnectivity(DM dm,PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes) 50 { 51 PetscErrorCode ierr; 52 PetscVTKInt *conn,*offsets; 53 PetscVTKType *types; 54 PetscInt dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,cMax,numLabelCells,hasLabel,c,v,countcell,countconn; 55 56 PetscFunctionBegin; 57 ierr = PetscMalloc3(piece->nconn,&conn,piece->ncells,&offsets,piece->ncells,&types);CHKERRQ(ierr); 58 59 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 60 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 61 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 62 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 63 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 64 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 65 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 66 ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr); 67 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 68 69 countcell = 0; 70 countconn = 0; 71 for (c = cStart; c < cEnd; ++c) { 72 PetscInt *closure = NULL; 73 PetscInt closureSize,nverts,celltype,startoffset,nC=0; 74 75 if (hasLabel) { 76 PetscInt value; 77 78 ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr); 79 if (value != 1) continue; 80 } 81 startoffset = countconn; 82 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 83 for (v = 0; v < closureSize*2; v += 2) { 84 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 85 conn[countconn++] = closure[v] - vStart; 86 ++nC; 87 } 88 } 89 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 90 ierr = DMPlexInvertCell(dim, nC, &conn[countconn-nC]);CHKERRQ(ierr); 91 92 offsets[countcell] = countconn; 93 94 nverts = countconn - startoffset; 95 ierr = DMPlexVTKGetCellType(dm,dim,nverts,&celltype);CHKERRQ(ierr); 96 97 types[countcell] = celltype; 98 countcell++; 99 } 100 if (countcell != piece->ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count"); 101 if (countconn != piece->nconn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count"); 102 *oconn = conn; 103 *ooffsets = offsets; 104 *otypes = types; 105 PetscFunctionReturn(0); 106 } 107 108 /* 109 Write all fields that have been provided to the viewer 110 Multi-block XML format with binary appended data. 111 */ 112 PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer) 113 { 114 MPI_Comm comm; 115 PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data; 116 PetscViewerVTKObjectLink link; 117 FILE *fp; 118 PetscMPIInt rank,size,tag; 119 PetscErrorCode ierr; 120 PetscInt dimEmbed,cellHeight,cStart,cEnd,vStart,vEnd,cMax,numLabelCells,hasLabel,c,v,r,i; 121 PieceInfo piece,*gpiece = NULL; 122 void *buffer = NULL; 123 124 PetscFunctionBegin; 125 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 126 #if defined(PETSC_USE_COMPLEX) 127 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Complex values not supported"); 128 #endif 129 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 130 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 131 tag = ((PetscObject)viewer)->tag; 132 133 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 134 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 135 #if defined(PETSC_WORDS_BIGENDIAN) 136 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr); 137 #else 138 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr); 139 #endif 140 ierr = PetscFPrintf(comm,fp," <UnstructuredGrid>\n");CHKERRQ(ierr); 141 142 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 143 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 144 ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr); 145 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 146 ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr); 147 if (cMax >= 0) cEnd = PetscMin(cEnd, cMax); 148 ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr); 149 150 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 151 piece.nvertices = vEnd - vStart; 152 piece.ncells = 0; 153 piece.nconn = 0; 154 for (c = cStart; c < cEnd; ++c) { 155 PetscInt *closure = NULL; 156 PetscInt closureSize; 157 158 if (hasLabel) { 159 PetscInt value; 160 161 ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr); 162 if (value != 1) continue; 163 } 164 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 165 for (v = 0; v < closureSize*2; v += 2) { 166 if ((closure[v] >= vStart) && (closure[v] < vEnd)) piece.nconn++; 167 } 168 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 169 piece.ncells++; 170 } 171 if (!rank) {ierr = PetscMalloc1(size,&gpiece);CHKERRQ(ierr);} 172 ierr = MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);CHKERRQ(ierr); 173 174 /* 175 * Write file header 176 */ 177 if (!rank) { 178 PetscInt boffset = 0; 179 180 for (r=0; r<size; r++) { 181 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);CHKERRQ(ierr); 182 /* Coordinate positions */ 183 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <Points>\n");CHKERRQ(ierr); 184 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 185 boffset += gpiece[r].nvertices*3*sizeof(PetscScalar) + sizeof(int); 186 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </Points>\n");CHKERRQ(ierr); 187 /* Cell connectivity */ 188 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <Cells>\n");CHKERRQ(ierr); 189 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 190 boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int); 191 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 192 boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int); 193 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 194 boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int); 195 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </Cells>\n");CHKERRQ(ierr); 196 197 /* 198 * Cell Data headers 199 */ 200 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <CellData>\n");CHKERRQ(ierr); 201 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 202 boffset += gpiece[r].ncells*sizeof(int) + sizeof(int); 203 /* all the vectors */ 204 for (link=vtk->link; link; link=link->next) { 205 Vec X = (Vec)link->vec; 206 PetscInt bs,nfields,field; 207 const char *vecname = ""; 208 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 209 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. */ 210 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 211 } 212 ierr = PetscSectionGetDof(dm->defaultSection,cStart,&bs);CHKERRQ(ierr); 213 ierr = PetscSectionGetNumFields(dm->defaultSection,&nfields);CHKERRQ(ierr); 214 for (field=0,i=0; field<(nfields?nfields:1); field++) { 215 PetscInt fbs,j; 216 PetscFV fv = NULL; 217 PetscObject f; 218 PetscClassId fClass; 219 const char *fieldname = NULL; 220 char buf[256]; 221 if (nfields) { /* We have user-defined fields/components */ 222 ierr = PetscSectionGetFieldDof(dm->defaultSection,cStart,field,&fbs);CHKERRQ(ierr); 223 ierr = PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);CHKERRQ(ierr); 224 } else fbs = bs; /* Say we have one field with 'bs' components */ 225 ierr = DMGetField(dm,field,&f);CHKERRQ(ierr); 226 ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr); 227 if (fClass == PETSCFV_CLASSID) { 228 fv = (PetscFV) f; 229 } 230 if (!fieldname) { 231 ierr = PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);CHKERRQ(ierr); 232 fieldname = buf; 233 } 234 for (j=0; j<fbs; j++) { 235 const char *compName = NULL; 236 if (fv) { 237 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr); 238 } 239 if (compName) { 240 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,compName,boffset);CHKERRQ(ierr); 241 } 242 else { 243 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);CHKERRQ(ierr); 244 } 245 boffset += gpiece[r].ncells*sizeof(PetscScalar) + sizeof(int); 246 i++; 247 } 248 } 249 if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs); 250 } 251 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n");CHKERRQ(ierr); 252 253 /* 254 * Point Data headers 255 */ 256 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n");CHKERRQ(ierr); 257 for (link=vtk->link; link; link=link->next) { 258 Vec X = (Vec)link->vec; 259 PetscInt bs,nfields,field; 260 const char *vecname = ""; 261 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 262 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. */ 263 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 264 } 265 ierr = PetscSectionGetDof(dm->defaultSection,vStart,&bs);CHKERRQ(ierr); 266 ierr = PetscSectionGetNumFields(dm->defaultSection,&nfields);CHKERRQ(ierr); 267 for (field=0,i=0; field<(nfields?nfields:1); field++) { 268 PetscInt fbs,j; 269 const char *fieldname = NULL; 270 char buf[256]; 271 if (nfields) { /* We have user-defined fields/components */ 272 ierr = PetscSectionGetFieldDof(dm->defaultSection,vStart,field,&fbs);CHKERRQ(ierr); 273 ierr = PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);CHKERRQ(ierr); 274 } else fbs = bs; /* Say we have one field with 'bs' components */ 275 if (!fieldname) { 276 ierr = PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);CHKERRQ(ierr); 277 fieldname = buf; 278 } 279 for (j=0; j<fbs; j++) { 280 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);CHKERRQ(ierr); 281 boffset += gpiece[r].nvertices*sizeof(PetscScalar) + sizeof(int); 282 } 283 } 284 } 285 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n");CHKERRQ(ierr); 286 287 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n");CHKERRQ(ierr); 288 } 289 } 290 291 ierr = PetscFPrintf(comm,fp," </UnstructuredGrid>\n");CHKERRQ(ierr); 292 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 293 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 294 295 if (!rank) { 296 PetscInt maxsize = 0; 297 for (r=0; r<size; r++) { 298 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscScalar))); 299 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*sizeof(PetscScalar))); 300 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt))); 301 } 302 ierr = PetscMalloc(maxsize,&buffer);CHKERRQ(ierr); 303 } 304 for (r=0; r<size; r++) { 305 if (r == rank) { 306 PetscInt nsend; 307 { /* Position */ 308 const PetscScalar *x; 309 PetscScalar *y = NULL; 310 Vec coords; 311 nsend = piece.nvertices*3; 312 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 313 ierr = VecGetArrayRead(coords,&x);CHKERRQ(ierr); 314 if (dimEmbed != 3) { 315 ierr = PetscMalloc1(piece.nvertices*3,&y);CHKERRQ(ierr); 316 for (i=0; i<piece.nvertices; i++) { 317 y[i*3+0] = x[i*dimEmbed+0]; 318 y[i*3+1] = (dimEmbed > 1) ? x[i*dimEmbed+1] : 0; 319 y[i*3+2] = 0; 320 } 321 } 322 ierr = TransferWrite(viewer,fp,r,0,y ? y : x,buffer,nsend,PETSC_SCALAR,tag);CHKERRQ(ierr); 323 ierr = PetscFree(y);CHKERRQ(ierr); 324 ierr = VecRestoreArrayRead(coords,&x);CHKERRQ(ierr); 325 } 326 { /* Connectivity, offsets, types */ 327 PetscVTKInt *connectivity = NULL, *offsets = NULL; 328 PetscVTKType *types = NULL; 329 ierr = DMPlexGetVTKConnectivity(dm,&piece,&connectivity,&offsets,&types);CHKERRQ(ierr); 330 ierr = TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,PETSC_INT32,tag);CHKERRQ(ierr); 331 ierr = TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,PETSC_INT32,tag);CHKERRQ(ierr); 332 ierr = TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,PETSC_UINT8,tag);CHKERRQ(ierr); 333 ierr = PetscFree3(connectivity,offsets,types);CHKERRQ(ierr); 334 } 335 { /* Owners (cell data) */ 336 PetscVTKInt *owners; 337 ierr = PetscMalloc1(piece.ncells,&owners);CHKERRQ(ierr); 338 for (i=0; i<piece.ncells; i++) owners[i] = rank; 339 ierr = TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,PETSC_INT32,tag);CHKERRQ(ierr); 340 ierr = PetscFree(owners);CHKERRQ(ierr); 341 } 342 /* Cell data */ 343 for (link=vtk->link; link; link=link->next) { 344 Vec X = (Vec)link->vec; 345 const PetscScalar *x; 346 PetscScalar *y; 347 PetscInt bs; 348 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 349 ierr = PetscSectionGetDof(dm->defaultSection,cStart,&bs);CHKERRQ(ierr); 350 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 351 ierr = PetscMalloc1(piece.ncells,&y);CHKERRQ(ierr); 352 for (i=0; i<bs; i++) { 353 PetscInt cnt; 354 for (c=cStart,cnt=0; c<cEnd; c++) { 355 PetscScalar *xpoint; 356 if (hasLabel) { /* Ignore some cells */ 357 PetscInt value; 358 ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr); 359 if (value != 1) continue; 360 } 361 ierr = DMPlexPointLocalRead(dm,c,x,&xpoint);CHKERRQ(ierr); 362 y[cnt++] = xpoint[i]; 363 } 364 if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 365 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,PETSC_SCALAR,tag);CHKERRQ(ierr); 366 } 367 ierr = PetscFree(y);CHKERRQ(ierr); 368 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 369 } 370 371 for (link=vtk->link; link; link=link->next) { 372 Vec X = (Vec)link->vec; 373 const PetscScalar *x; 374 PetscScalar *y; 375 PetscInt bs; 376 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 377 ierr = PetscSectionGetDof(dm->defaultSection,vStart,&bs);CHKERRQ(ierr); 378 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 379 ierr = PetscMalloc1(piece.nvertices,&y);CHKERRQ(ierr); 380 for (i=0; i<bs; i++) { 381 PetscInt cnt; 382 for (v=vStart,cnt=0; v<vEnd; v++) { 383 PetscScalar *xpoint; 384 ierr = DMPlexPointLocalRead(dm,v,x,&xpoint);CHKERRQ(ierr); 385 y[cnt++] = xpoint[i]; 386 } 387 if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 388 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,PETSC_SCALAR,tag);CHKERRQ(ierr); 389 } 390 ierr = PetscFree(y);CHKERRQ(ierr); 391 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 392 } 393 } else if (!rank) { 394 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,PETSC_SCALAR,tag);CHKERRQ(ierr); /* positions */ 395 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,PETSC_INT32,tag);CHKERRQ(ierr); /* connectivity */ 396 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag);CHKERRQ(ierr); /* offsets */ 397 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_UINT8,tag);CHKERRQ(ierr); /* types */ 398 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag);CHKERRQ(ierr); /* owner rank (cells) */ 399 /* all cell data */ 400 for (link=vtk->link; link; link=link->next) { 401 PetscInt bs; 402 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 403 ierr = PetscSectionGetDof(dm->defaultSection,cStart,&bs);CHKERRQ(ierr); 404 for (i=0; i<bs; i++) { 405 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_SCALAR,tag);CHKERRQ(ierr); 406 } 407 } 408 /* all point data */ 409 for (link=vtk->link; link; link=link->next) { 410 PetscInt bs; 411 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 412 ierr = PetscSectionGetDof(dm->defaultSection,vStart,&bs);CHKERRQ(ierr); 413 for (i=0; i<bs; i++) { 414 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,PETSC_SCALAR,tag);CHKERRQ(ierr); 415 } 416 } 417 } 418 } 419 ierr = PetscFree(gpiece);CHKERRQ(ierr); 420 ierr = PetscFree(buffer);CHKERRQ(ierr); 421 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 422 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 423 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 424 PetscFunctionReturn(0); 425 } 426