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