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 const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian"; 133 134 PetscFunctionBegin; 135 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 136 #if defined(PETSC_USE_COMPLEX) 137 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Complex values not supported"); 138 #endif 139 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 140 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 141 tag = ((PetscObject)viewer)->tag; 142 143 ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr); 144 ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr); 145 ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);CHKERRQ(ierr); 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 = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr); 153 ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr); 154 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 155 156 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 157 piece.nvertices = 0; 158 piece.ncells = 0; 159 piece.nconn = 0; 160 if (!localized) piece.nvertices = vEnd - vStart; 161 for (c = cStart; c < cEnd; ++c) { 162 PetscInt dof = 0; 163 if (hasLabel) { 164 PetscInt value; 165 166 ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr); 167 if (value != 1) continue; 168 } 169 if (localized) { 170 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 171 } 172 if (!dof) { 173 PetscInt *closure = NULL; 174 PetscInt closureSize; 175 176 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 177 for (v = 0; v < closureSize*2; v += 2) { 178 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 179 piece.nconn++; 180 if (localized) piece.nvertices++; 181 } 182 } 183 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 184 } else { 185 piece.nvertices += dof/dimEmbed; 186 piece.nconn += dof/dimEmbed; 187 } 188 piece.ncells++; 189 } 190 if (!rank) {ierr = PetscMalloc1(size,&gpiece);CHKERRQ(ierr);} 191 ierr = MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);CHKERRQ(ierr); 192 193 /* 194 * Write file header 195 */ 196 if (!rank) { 197 PetscInt boffset = 0; 198 199 for (r=0; r<size; r++) { 200 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);CHKERRQ(ierr); 201 /* Coordinate positions */ 202 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <Points>\n");CHKERRQ(ierr); 203 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr); 204 boffset += gpiece[r].nvertices*3*sizeof(PetscScalar) + sizeof(int); 205 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </Points>\n");CHKERRQ(ierr); 206 /* Cell connectivity */ 207 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <Cells>\n");CHKERRQ(ierr); 208 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 209 boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int); 210 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 211 boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int); 212 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 213 boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int); 214 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </Cells>\n");CHKERRQ(ierr); 215 216 /* 217 * Cell Data headers 218 */ 219 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <CellData>\n");CHKERRQ(ierr); 220 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr); 221 boffset += gpiece[r].ncells*sizeof(int) + sizeof(int); 222 /* all the vectors */ 223 for (link=vtk->link; link; link=link->next) { 224 Vec X = (Vec)link->vec; 225 DM dmX = NULL; 226 PetscInt bs,nfields,field; 227 const char *vecname = ""; 228 PetscSection section; 229 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 230 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. */ 231 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 232 } 233 ierr = VecGetDM(X, &dmX);CHKERRQ(ierr); 234 if (!dmX) dmX = dm; 235 ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);CHKERRQ(ierr); 236 if (!section) {ierr = DMGetLocalSection(dmX, §ion);CHKERRQ(ierr);} 237 ierr = PetscSectionGetDof(section,cStart,&bs);CHKERRQ(ierr); 238 ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr); 239 field = 0; 240 if (link->field >= 0) { 241 field = link->field; 242 nfields = field + 1; 243 } 244 for (i=0; field<(nfields?nfields:1); field++) { 245 PetscInt fbs,j; 246 PetscFV fv = NULL; 247 PetscObject f; 248 PetscClassId fClass; 249 const char *fieldname = NULL; 250 char buf[256]; 251 PetscBool vector; 252 if (nfields) { /* We have user-defined fields/components */ 253 ierr = PetscSectionGetFieldDof(section,cStart,field,&fbs);CHKERRQ(ierr); 254 ierr = PetscSectionGetFieldName(section,field,&fieldname);CHKERRQ(ierr); 255 } else fbs = bs; /* Say we have one field with 'bs' components */ 256 ierr = DMGetField(dmX,field,NULL,&f);CHKERRQ(ierr); 257 ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr); 258 if (fClass == PETSCFV_CLASSID) { 259 fv = (PetscFV) f; 260 } 261 if (nfields && !fieldname) { 262 ierr = PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);CHKERRQ(ierr); 263 fieldname = buf; 264 } 265 vector = PETSC_FALSE; 266 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 267 vector = PETSC_TRUE; 268 if (fbs > 3) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Cell vector fields can have at most 3 components, %D given\n", fbs); 269 for (j = 0; j < fbs; j++) { 270 const char *compName = NULL; 271 if (fv) { 272 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr); 273 if (compName) break; 274 } 275 } 276 if (j < fbs) vector = PETSC_FALSE; 277 } 278 if (vector) { 279 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 280 boffset += gpiece[r].ncells*3*sizeof(PetscScalar) + sizeof(int); 281 i+=fbs; 282 } else { 283 for (j=0; j<fbs; j++) { 284 const char *compName = NULL; 285 if (fv) { 286 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr); 287 } 288 if (compName) { 289 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); 290 } 291 else if (fbs > 1) { 292 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); 293 } else { 294 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 295 } 296 boffset += gpiece[r].ncells*sizeof(PetscScalar) + sizeof(int); 297 i++; 298 } 299 } 300 } 301 if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs); 302 } 303 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n");CHKERRQ(ierr); 304 305 /* 306 * Point Data headers 307 */ 308 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n");CHKERRQ(ierr); 309 for (link=vtk->link; link; link=link->next) { 310 Vec X = (Vec)link->vec; 311 DM dmX; 312 PetscInt bs,nfields,field; 313 const char *vecname = ""; 314 PetscSection section; 315 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 316 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. */ 317 ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr); 318 } 319 ierr = VecGetDM(X, &dmX);CHKERRQ(ierr); 320 if (!dmX) dmX = dm; 321 ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);CHKERRQ(ierr); 322 if (!section) {ierr = DMGetLocalSection(dmX, §ion);CHKERRQ(ierr);} 323 ierr = PetscSectionGetDof(section,vStart,&bs);CHKERRQ(ierr); 324 ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr); 325 field = 0; 326 if (link->field >= 0) { 327 field = link->field; 328 nfields = field + 1; 329 } 330 for (i=0; field<(nfields?nfields:1); field++) { 331 PetscInt fbs,j; 332 const char *fieldname = NULL; 333 char buf[256]; 334 if (nfields) { /* We have user-defined fields/components */ 335 ierr = PetscSectionGetFieldDof(section,vStart,field,&fbs);CHKERRQ(ierr); 336 ierr = PetscSectionGetFieldName(section,field,&fieldname);CHKERRQ(ierr); 337 } else fbs = bs; /* Say we have one field with 'bs' components */ 338 if (nfields && !fieldname) { 339 ierr = PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);CHKERRQ(ierr); 340 fieldname = buf; 341 } 342 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 343 if (fbs > 3) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Point vector fields can have at most 3 components, %D given\n", fbs); 344 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 345 boffset += gpiece[r].nvertices*3*sizeof(PetscScalar) + sizeof(int); 346 } else { 347 for (j=0; j<fbs; j++) { 348 if (fbs > 1) { 349 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); 350 } else { 351 ierr = PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr); 352 } 353 boffset += gpiece[r].nvertices*sizeof(PetscScalar) + sizeof(int); 354 } 355 } 356 } 357 } 358 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n");CHKERRQ(ierr); 359 ierr = PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n");CHKERRQ(ierr); 360 } 361 } 362 363 ierr = PetscFPrintf(comm,fp," </UnstructuredGrid>\n");CHKERRQ(ierr); 364 ierr = PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr); 365 ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr); 366 367 if (!rank) { 368 PetscInt maxsize = 0; 369 for (r=0; r<size; r++) { 370 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscScalar))); 371 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*3*sizeof(PetscScalar))); 372 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt))); 373 } 374 ierr = PetscMalloc(maxsize,&buffer);CHKERRQ(ierr); 375 } 376 for (r=0; r<size; r++) { 377 if (r == rank) { 378 PetscInt nsend; 379 { /* Position */ 380 const PetscScalar *x; 381 PetscScalar *y = NULL; 382 Vec coords; 383 384 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 385 ierr = VecGetArrayRead(coords,&x);CHKERRQ(ierr); 386 if (dimEmbed != 3 || localized) { 387 ierr = PetscMalloc1(piece.nvertices*3,&y);CHKERRQ(ierr); 388 if (localized) { 389 PetscInt cnt; 390 for (c=cStart,cnt=0; c<cEnd; c++) { 391 PetscInt off, dof; 392 393 ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr); 394 if (!dof) { 395 PetscInt *closure = NULL; 396 PetscInt closureSize; 397 398 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 399 for (v = 0; v < closureSize*2; v += 2) { 400 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 401 ierr = PetscSectionGetOffset(coordSection, closure[v], &off);CHKERRQ(ierr); 402 if (dimEmbed != 3) { 403 y[cnt*3+0] = x[off+0]; 404 y[cnt*3+1] = (dimEmbed > 1) ? x[off+1] : 0.0; 405 y[cnt*3+2] = 0.0; 406 } else { 407 y[cnt*3+0] = x[off+0]; 408 y[cnt*3+1] = x[off+1]; 409 y[cnt*3+2] = x[off+2]; 410 } 411 cnt++; 412 } 413 } 414 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 415 } else { 416 ierr = PetscSectionGetOffset(coordSection, c, &off);CHKERRQ(ierr); 417 if (dimEmbed != 3) { 418 for (i=0; i<dof/dimEmbed; i++) { 419 y[cnt*3+0] = x[off + i*dimEmbed + 0]; 420 y[cnt*3+1] = (dimEmbed > 1) ? x[off + i*dimEmbed + 1] : 0.0; 421 y[cnt*3+2] = 0.0; 422 cnt++; 423 } 424 } else { 425 for (i=0; i<dof; i ++) { 426 y[cnt*3+i] = x[off + i]; 427 } 428 cnt += dof/dimEmbed; 429 } 430 } 431 } 432 if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 433 } else { 434 for (i=0; i<piece.nvertices; i++) { 435 y[i*3+0] = x[i*dimEmbed+0]; 436 y[i*3+1] = (dimEmbed > 1) ? x[i*dimEmbed+1] : 0; 437 y[i*3+2] = 0.0; 438 } 439 } 440 } 441 nsend = piece.nvertices*3; 442 ierr = TransferWrite(viewer,fp,r,0,y ? y : x,buffer,nsend,MPIU_SCALAR,tag);CHKERRQ(ierr); 443 ierr = PetscFree(y);CHKERRQ(ierr); 444 ierr = VecRestoreArrayRead(coords,&x);CHKERRQ(ierr); 445 } 446 { /* Connectivity, offsets, types */ 447 PetscVTKInt *connectivity = NULL, *offsets = NULL; 448 PetscVTKType *types = NULL; 449 ierr = DMPlexGetVTKConnectivity(dm,localized,&piece,&connectivity,&offsets,&types);CHKERRQ(ierr); 450 ierr = TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,MPI_INT,tag);CHKERRQ(ierr); 451 ierr = TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,MPI_INT,tag);CHKERRQ(ierr); 452 ierr = TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,MPI_CHAR,tag);CHKERRQ(ierr); 453 ierr = PetscFree3(connectivity,offsets,types);CHKERRQ(ierr); 454 } 455 { /* Owners (cell data) */ 456 PetscVTKInt *owners; 457 ierr = PetscMalloc1(piece.ncells,&owners);CHKERRQ(ierr); 458 for (i=0; i<piece.ncells; i++) owners[i] = rank; 459 ierr = TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,MPI_INT,tag);CHKERRQ(ierr); 460 ierr = PetscFree(owners);CHKERRQ(ierr); 461 } 462 /* Cell data */ 463 for (link=vtk->link; link; link=link->next) { 464 Vec X = (Vec)link->vec; 465 DM dmX; 466 const PetscScalar *x; 467 PetscScalar *y; 468 PetscInt bs, nfields, field; 469 PetscSection section = NULL; 470 471 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 472 ierr = VecGetDM(X, &dmX);CHKERRQ(ierr); 473 if (!dmX) dmX = dm; 474 ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);CHKERRQ(ierr); 475 if (!section) {ierr = DMGetLocalSection(dmX, §ion);CHKERRQ(ierr);} 476 ierr = PetscSectionGetDof(section,cStart,&bs);CHKERRQ(ierr); 477 ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr); 478 field = 0; 479 if (link->field >= 0) { 480 field = link->field; 481 nfields = field + 1; 482 } 483 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 484 ierr = PetscMalloc1(piece.ncells*3,&y);CHKERRQ(ierr); 485 for (i=0; field<(nfields?nfields:1); field++) { 486 PetscInt fbs,j; 487 PetscFV fv = NULL; 488 PetscObject f; 489 PetscClassId fClass; 490 PetscBool vector; 491 if (nfields) { /* We have user-defined fields/components */ 492 ierr = PetscSectionGetFieldDof(section,cStart,field,&fbs);CHKERRQ(ierr); 493 } else fbs = bs; /* Say we have one field with 'bs' components */ 494 ierr = DMGetField(dmX,field,NULL,&f);CHKERRQ(ierr); 495 ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr); 496 if (fClass == PETSCFV_CLASSID) { 497 fv = (PetscFV) f; 498 } 499 vector = PETSC_FALSE; 500 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 501 vector = PETSC_TRUE; 502 for (j = 0; j < fbs; j++) { 503 const char *compName = NULL; 504 if (fv) { 505 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr); 506 if (compName) break; 507 } 508 } 509 if (j < fbs) vector = PETSC_FALSE; 510 } 511 if (vector) { 512 PetscInt cnt; 513 for (c=cStart,cnt=0; c<cEnd; c++) { 514 const PetscScalar *xpoint; 515 PetscInt off, j; 516 517 if (hasLabel) { /* Ignore some cells */ 518 PetscInt value; 519 ierr = DMGetLabelValue(dmX, "vtk", c, &value);CHKERRQ(ierr); 520 if (value != 1) continue; 521 } 522 if (nfields) { 523 ierr = PetscSectionGetFieldOffset(section, c, field, &off);CHKERRQ(ierr); 524 } else { 525 ierr = PetscSectionGetOffset(section, c, &off);CHKERRQ(ierr); 526 } 527 xpoint = &x[off]; 528 for (j = 0; j < fbs; j++) { 529 y[cnt++] = xpoint[j]; 530 } 531 for (; j < 3; j++) y[cnt++] = 0.; 532 } 533 if (cnt != piece.ncells*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 534 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells*3,MPIU_SCALAR,tag);CHKERRQ(ierr); 535 } else { 536 for (i=0; i<fbs; i++) { 537 PetscInt cnt; 538 for (c=cStart,cnt=0; c<cEnd; c++) { 539 const PetscScalar *xpoint; 540 PetscInt off; 541 542 if (hasLabel) { /* Ignore some cells */ 543 PetscInt value; 544 ierr = DMGetLabelValue(dmX, "vtk", c, &value);CHKERRQ(ierr); 545 if (value != 1) continue; 546 } 547 if (nfields) { 548 ierr = PetscSectionGetFieldOffset(section, c, field, &off);CHKERRQ(ierr); 549 } else { 550 ierr = PetscSectionGetOffset(section, c, &off);CHKERRQ(ierr); 551 } 552 xpoint = &x[off]; 553 y[cnt++] = xpoint[i]; 554 } 555 if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 556 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,MPIU_SCALAR,tag);CHKERRQ(ierr); 557 } 558 } 559 } 560 ierr = PetscFree(y);CHKERRQ(ierr); 561 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 562 } 563 /* point data */ 564 for (link=vtk->link; link; link=link->next) { 565 Vec X = (Vec)link->vec; 566 DM dmX; 567 const PetscScalar *x; 568 PetscScalar *y; 569 PetscInt bs, nfields, field; 570 PetscSection section = NULL; 571 572 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 573 ierr = VecGetDM(X, &dmX);CHKERRQ(ierr); 574 if (!dmX) dmX = dm; 575 ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);CHKERRQ(ierr); 576 if (!section) {ierr = DMGetLocalSection(dmX, §ion);CHKERRQ(ierr);} 577 ierr = PetscSectionGetDof(section,vStart,&bs);CHKERRQ(ierr); 578 ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr); 579 field = 0; 580 if (link->field >= 0) { 581 field = link->field; 582 nfields = field + 1; 583 } 584 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 585 ierr = PetscMalloc1(piece.nvertices*3,&y);CHKERRQ(ierr); 586 for (i=0; field<(nfields?nfields:1); field++) { 587 PetscInt fbs,j; 588 if (nfields) { /* We have user-defined fields/components */ 589 ierr = PetscSectionGetFieldDof(section,vStart,field,&fbs);CHKERRQ(ierr); 590 } else fbs = bs; /* Say we have one field with 'bs' components */ 591 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 592 PetscInt cnt; 593 if (!localized) { 594 for (v=vStart,cnt=0; v<vEnd; v++) { 595 PetscInt off; 596 const PetscScalar *xpoint; 597 598 if (nfields) { 599 ierr = PetscSectionGetFieldOffset(section,v,field,&off);CHKERRQ(ierr); 600 } else { 601 ierr = PetscSectionGetOffset(section,v,&off);CHKERRQ(ierr); 602 } 603 xpoint = &x[off]; 604 for (j = 0; j < fbs; j++) { 605 y[cnt++] = xpoint[j]; 606 } 607 for (; j < 3; j++) y[cnt++] = 0.; 608 } 609 } else { 610 for (c=cStart,cnt=0; c<cEnd; c++) { 611 PetscInt *closure = NULL; 612 PetscInt closureSize, off; 613 614 ierr = DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 615 for (v = 0, off = 0; v < closureSize*2; v += 2) { 616 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 617 PetscInt voff; 618 const PetscScalar *xpoint; 619 620 if (nfields) { 621 ierr = PetscSectionGetFieldOffset(section,closure[v],field,&voff);CHKERRQ(ierr); 622 } else { 623 ierr = PetscSectionGetOffset(section,closure[v],&voff);CHKERRQ(ierr); 624 } 625 xpoint = &x[voff]; 626 for (j = 0; j < fbs; j++) { 627 y[cnt + off++] = xpoint[i]; 628 } 629 for (; j < 3; j++) y[cnt + off++] = 0.; 630 } 631 } 632 cnt += off; 633 ierr = DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 634 } 635 } 636 if (cnt != piece.nvertices*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 637 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices*3,MPIU_SCALAR,tag);CHKERRQ(ierr); 638 } else { 639 for (i=0; i<fbs; i++) { 640 PetscInt cnt; 641 if (!localized) { 642 for (v=vStart,cnt=0; v<vEnd; v++) { 643 PetscInt off; 644 const PetscScalar *xpoint; 645 646 if (nfields) { 647 ierr = PetscSectionGetFieldOffset(section,v,field,&off);CHKERRQ(ierr); 648 } else { 649 ierr = PetscSectionGetOffset(section,v,&off);CHKERRQ(ierr); 650 } 651 xpoint = &x[off]; 652 y[cnt++] = xpoint[i]; 653 } 654 } else { 655 for (c=cStart,cnt=0; c<cEnd; c++) { 656 PetscInt *closure = NULL; 657 PetscInt closureSize, off; 658 659 ierr = DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 660 for (v = 0, off = 0; v < closureSize*2; v += 2) { 661 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 662 PetscInt voff; 663 const PetscScalar *xpoint; 664 665 if (nfields) { 666 ierr = PetscSectionGetFieldOffset(section,closure[v],field,&voff);CHKERRQ(ierr); 667 } else { 668 ierr = PetscSectionGetOffset(section,closure[v],&voff);CHKERRQ(ierr); 669 } 670 xpoint = &x[voff]; 671 y[cnt + off++] = xpoint[i]; 672 } 673 } 674 cnt += off; 675 ierr = DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 676 } 677 } 678 if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 679 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,MPIU_SCALAR,tag);CHKERRQ(ierr); 680 } 681 } 682 } 683 ierr = PetscFree(y);CHKERRQ(ierr); 684 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 685 } 686 } else if (!rank) { 687 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_SCALAR,tag);CHKERRQ(ierr); /* positions */ 688 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,MPI_INT,tag);CHKERRQ(ierr); /* connectivity */ 689 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag);CHKERRQ(ierr); /* offsets */ 690 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_CHAR,tag);CHKERRQ(ierr); /* types */ 691 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag);CHKERRQ(ierr); /* owner rank (cells) */ 692 /* all cell data */ 693 for (link=vtk->link; link; link=link->next) { 694 Vec X = (Vec)link->vec; 695 PetscInt bs, nfields, field; 696 DM dmX; 697 PetscSection section = NULL; 698 699 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 700 ierr = VecGetDM(X, &dmX);CHKERRQ(ierr); 701 if (!dmX) dmX = dm; 702 ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);CHKERRQ(ierr); 703 if (!section) {ierr = DMGetLocalSection(dmX, §ion);CHKERRQ(ierr);} 704 ierr = PetscSectionGetDof(section,cStart,&bs);CHKERRQ(ierr); 705 ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr); 706 field = 0; 707 if (link->field >= 0) { 708 field = link->field; 709 nfields = field + 1; 710 } 711 for (i=0; field<(nfields?nfields:1); field++) { 712 PetscInt fbs,j; 713 PetscFV fv = NULL; 714 PetscObject f; 715 PetscClassId fClass; 716 PetscBool vector; 717 if (nfields) { /* We have user-defined fields/components */ 718 ierr = PetscSectionGetFieldDof(section,cStart,field,&fbs);CHKERRQ(ierr); 719 } else fbs = bs; /* Say we have one field with 'bs' components */ 720 ierr = DMGetField(dmX,field,NULL,&f);CHKERRQ(ierr); 721 ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr); 722 if (fClass == PETSCFV_CLASSID) { 723 fv = (PetscFV) f; 724 } 725 vector = PETSC_FALSE; 726 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 727 vector = PETSC_TRUE; 728 for (j = 0; j < fbs; j++) { 729 const char *compName = NULL; 730 if (fv) { 731 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr); 732 if (compName) break; 733 } 734 } 735 if (j < fbs) vector = PETSC_FALSE; 736 } 737 if (vector) { 738 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells*3,MPIU_SCALAR,tag);CHKERRQ(ierr); 739 } else { 740 for (i=0; i<fbs; i++) { 741 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPIU_SCALAR,tag);CHKERRQ(ierr); 742 } 743 } 744 } 745 } 746 /* all point data */ 747 for (link=vtk->link; link; link=link->next) { 748 Vec X = (Vec)link->vec; 749 DM dmX; 750 PetscInt bs, nfields, field; 751 PetscSection section = NULL; 752 753 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 754 ierr = VecGetDM(X, &dmX);CHKERRQ(ierr); 755 if (!dmX) dmX = dm; 756 ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);CHKERRQ(ierr); 757 if (!section) {ierr = DMGetLocalSection(dmX, §ion);CHKERRQ(ierr);} 758 ierr = PetscSectionGetDof(section,vStart,&bs);CHKERRQ(ierr); 759 ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr); 760 field = 0; 761 if (link->field >= 0) { 762 field = link->field; 763 nfields = field + 1; 764 } 765 for (i=0; field<(nfields?nfields:1); field++) { 766 PetscInt fbs; 767 if (nfields) { /* We have user-defined fields/components */ 768 ierr = PetscSectionGetFieldDof(section,vStart,field,&fbs);CHKERRQ(ierr); 769 } else fbs = bs; /* Say we have one field with 'bs' components */ 770 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 771 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_SCALAR,tag);CHKERRQ(ierr); 772 } else { 773 for (i=0; i<fbs; i++) { 774 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,MPIU_SCALAR,tag);CHKERRQ(ierr); 775 } 776 } 777 } 778 } 779 } 780 } 781 ierr = PetscFree(gpiece);CHKERRQ(ierr); 782 ierr = PetscFree(buffer);CHKERRQ(ierr); 783 ierr = PetscFPrintf(comm,fp,"\n </AppendedData>\n");CHKERRQ(ierr); 784 ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr); 785 ierr = PetscFClose(comm,fp);CHKERRQ(ierr); 786 PetscFunctionReturn(0); 787 } 788