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