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