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