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) || defined(PETSC_USE_REAL___FP16) 11 /* output in float if single or half precision in memory */ 12 static const char precision[] = "Float32"; 13 typedef float PetscVTUReal; 14 #define MPIU_VTUREAL MPI_FLOAT 15 #elif defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 16 /* output in double if double or quad precision in memory */ 17 static const char precision[] = "Float64"; 18 typedef double PetscVTUReal; 19 #define MPIU_VTUREAL MPI_DOUBLE 20 #else 21 static const char precision[] = "UnknownPrecision"; 22 typedef PetscReal PetscVTUReal; 23 #define MPIU_VTUREAL MPIU_REAL 24 #endif 25 26 static PetscErrorCode TransferWrite(MPI_Comm comm, PetscViewer viewer,FILE *fp,PetscMPIInt srank,PetscMPIInt root,const void *send,void *recv,PetscMPIInt count,MPI_Datatype mpidatatype,PetscMPIInt tag) 27 { 28 PetscMPIInt rank; 29 30 PetscFunctionBegin; 31 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 32 if (rank == srank && rank != root) { 33 PetscCallMPI(MPI_Send((void*)send,count,mpidatatype,root,tag,comm)); 34 } else if (rank == root) { 35 const void *buffer; 36 if (root == srank) { /* self */ 37 buffer = send; 38 } else { 39 MPI_Status status; 40 PetscMPIInt nrecv; 41 PetscCallMPI(MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status)); 42 PetscCallMPI(MPI_Get_count(&status,mpidatatype,&nrecv)); 43 PetscCheck(count == nrecv,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch"); 44 buffer = recv; 45 } 46 PetscCall(PetscViewerVTKFWrite(viewer,fp,buffer,count,mpidatatype)); 47 } 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes) 52 { 53 PetscSection coordSection, cellCoordSection; 54 PetscVTKInt *conn,*offsets; 55 PetscVTKType *types; 56 PetscInt dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,numLabelCells,hasLabel,c,v,countcell,countconn; 57 58 PetscFunctionBegin; 59 PetscCall(PetscMalloc3(piece->nconn,&conn,piece->ncells,&offsets,piece->ncells,&types)); 60 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 61 PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection)); 62 PetscCall(DMGetDimension(dm,&dim)); 63 PetscCall(DMPlexGetChart(dm,&pStart,&pEnd)); 64 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 65 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 66 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 67 PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells)); 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 PetscCall(DMGetLabelValue(dm, "vtk", c, &value)); 79 if (value != 1) continue; 80 } 81 startoffset = countconn; 82 if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 83 if (!dof) { 84 PetscInt *closure = NULL; 85 PetscInt closureSize; 86 87 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 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 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 96 } else { 97 for (nC = 0; nC < dof/dim; nC++) conn[countconn++] = startoffset + nC; 98 } 99 100 { 101 PetscInt n = PetscMin(nC, 8), s = countconn - nC, i, cone[8]; 102 for (i = 0; i < n; ++i) cone[i] = conn[s+i]; 103 PetscCall(DMPlexReorderCell(dm, c, cone)); 104 for (i = 0; i < n; ++i) conn[s+i] = (int)cone[i]; 105 } 106 107 offsets[countcell] = countconn; 108 109 nverts = countconn - startoffset; 110 PetscCall(DMPlexVTKGetCellType_Internal(dm,dim,nverts,&celltype)); 111 112 types[countcell] = celltype; 113 countcell++; 114 } 115 PetscCheck(countcell == piece->ncells,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count"); 116 PetscCheck(countconn == piece->nconn,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count"); 117 *oconn = conn; 118 *ooffsets = offsets; 119 *otypes = types; 120 PetscFunctionReturn(0); 121 } 122 123 /* 124 Write all fields that have been provided to the viewer 125 Multi-block XML format with binary appended data. 126 */ 127 PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer) 128 { 129 MPI_Comm comm; 130 PetscSection coordSection, cellCoordSection; 131 PetscViewer_VTK *vtk = (PetscViewer_VTK*)viewer->data; 132 PetscViewerVTKObjectLink link; 133 FILE *fp; 134 PetscMPIInt rank,size,tag; 135 PetscInt dimEmbed,cellHeight,cStart,cEnd,vStart,vEnd,numLabelCells,hasLabel,c,v,r,i; 136 PetscBool localized; 137 PieceInfo piece,*gpiece = NULL; 138 void *buffer = NULL; 139 const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian"; 140 PetscInt loops_per_scalar; 141 142 PetscFunctionBegin; 143 PetscCall(PetscObjectGetComm((PetscObject)dm,&comm)); 144 #if defined(PETSC_USE_COMPLEX) 145 loops_per_scalar = 2; 146 #else 147 loops_per_scalar = 1; 148 #endif 149 PetscCallMPI(MPI_Comm_size(comm,&size)); 150 PetscCallMPI(MPI_Comm_rank(comm,&rank)); 151 PetscCall(PetscCommGetNewTag(comm,&tag)); 152 153 PetscCall(PetscFOpen(comm,vtk->filename,"wb",&fp)); 154 PetscCall(PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n")); 155 PetscCall(PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order)); 156 PetscCall(PetscFPrintf(comm,fp," <UnstructuredGrid>\n")); 157 158 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 159 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight)); 160 PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd)); 161 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 162 PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells)); 163 PetscCall(DMGetCoordinatesLocalized(dm, &localized)); 164 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 165 PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection)); 166 167 hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE; 168 piece.nvertices = 0; 169 piece.ncells = 0; 170 piece.nconn = 0; 171 if (!localized) piece.nvertices = vEnd - vStart; 172 for (c = cStart; c < cEnd; ++c) { 173 PetscInt dof = 0; 174 if (hasLabel) { 175 PetscInt value; 176 177 PetscCall(DMGetLabelValue(dm, "vtk", c, &value)); 178 if (value != 1) continue; 179 } 180 if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 181 if (!dof) { 182 PetscInt *closure = NULL; 183 PetscInt closureSize; 184 185 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 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 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 193 } else { 194 piece.nvertices += dof/dimEmbed; 195 piece.nconn += dof/dimEmbed; 196 } 197 piece.ncells++; 198 } 199 if (rank == 0) PetscCall(PetscMalloc1(size,&gpiece)); 200 PetscCallMPI(MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm)); 201 202 /* 203 * Write file header 204 */ 205 if (rank == 0) { 206 PetscInt boffset = 0; 207 208 for (r=0; r<size; r++) { 209 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n",gpiece[r].nvertices,gpiece[r].ncells)); 210 /* Coordinate positions */ 211 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <Points>\n")); 212 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,boffset)); 213 boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int); 214 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," </Points>\n")); 215 /* Cell connectivity */ 216 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <Cells>\n")); 217 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",boffset)); 218 boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int); 219 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",boffset)); 220 boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int); 221 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",boffset)); 222 boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int); 223 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," </Cells>\n")); 224 225 /* 226 * Cell Data headers 227 */ 228 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <CellData>\n")); 229 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",boffset)); 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 DM dmX = NULL; 235 PetscInt bs = 1,nfields,field; 236 const char *vecname = ""; 237 PetscSection section; 238 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 239 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. */ 240 PetscCall(PetscObjectGetName((PetscObject)X,&vecname)); 241 } 242 PetscCall(VecGetDM(X, &dmX)); 243 if (!dmX) dmX = dm; 244 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject*) §ion)); 245 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 246 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section,cStart,&bs)); 247 PetscCall(PetscSectionGetNumFields(section,&nfields)); 248 field = 0; 249 if (link->field >= 0) { 250 field = link->field; 251 nfields = field + 1; 252 } 253 for (i=0; field<(nfields?nfields:1); field++) { 254 PetscInt fbs,j; 255 PetscFV fv = NULL; 256 PetscObject f; 257 PetscClassId fClass; 258 const char *fieldname = NULL; 259 char buf[256]; 260 PetscBool vector; 261 if (nfields) { /* We have user-defined fields/components */ 262 PetscCall(PetscSectionGetFieldDof(section,cStart,field,&fbs)); 263 PetscCall(PetscSectionGetFieldName(section,field,&fieldname)); 264 } else fbs = bs; /* Say we have one field with 'bs' components */ 265 PetscCall(DMGetField(dmX,field,NULL,&f)); 266 PetscCall(PetscObjectGetClassId(f,&fClass)); 267 if (fClass == PETSCFV_CLASSID) { 268 fv = (PetscFV) f; 269 } 270 if (nfields && !fieldname) { 271 PetscCall(PetscSNPrintf(buf,sizeof(buf),"CellField%" PetscInt_FMT,field)); 272 fieldname = buf; 273 } 274 vector = PETSC_FALSE; 275 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 276 vector = PETSC_TRUE; 277 PetscCheck(fbs <= 3,PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Cell vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs); 278 for (j = 0; j < fbs; j++) { 279 const char *compName = NULL; 280 if (fv) { 281 PetscCall(PetscFVGetComponentName(fv,j,&compName)); 282 if (compName) break; 283 } 284 } 285 if (j < fbs) vector = PETSC_FALSE; 286 } 287 if (vector) { 288 #if defined(PETSC_USE_COMPLEX) 289 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,vecname,fieldname,boffset)); 290 boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int); 291 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,vecname,fieldname,boffset)); 292 boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int); 293 #else 294 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,vecname,fieldname,boffset)); 295 boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int); 296 #endif 297 i+=fbs; 298 } else { 299 for (j=0; j<fbs; j++) { 300 const char *compName = NULL; 301 char finalname[256]; 302 if (fv) { 303 PetscCall(PetscFVGetComponentName(fv,j,&compName)); 304 } 305 if (compName) { 306 PetscCall(PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName)); 307 } 308 else if (fbs > 1) { 309 PetscCall(PetscSNPrintf(finalname,255,"%s%s.%" PetscInt_FMT,vecname,fieldname,j)); 310 } else { 311 PetscCall(PetscSNPrintf(finalname,255,"%s%s",vecname,fieldname)); 312 } 313 #if defined(PETSC_USE_COMPLEX) 314 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,finalname,boffset)); 315 boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int); 316 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,finalname,boffset)); 317 boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int); 318 #else 319 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,finalname,boffset)); 320 boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int); 321 #endif 322 i++; 323 } 324 } 325 } 326 //PetscCheck(i == bs,comm,PETSC_ERR_PLIB,"Total number of field components %" PetscInt_FMT " != block size %" PetscInt_FMT,i,bs); 327 } 328 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," </CellData>\n")); 329 330 /* 331 * Point Data headers 332 */ 333 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," <PointData>\n")); 334 for (link=vtk->link; link; link=link->next) { 335 Vec X = (Vec)link->vec; 336 DM dmX; 337 PetscInt bs = 1,nfields,field; 338 const char *vecname = ""; 339 PetscSection section; 340 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 341 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. */ 342 PetscCall(PetscObjectGetName((PetscObject)X,&vecname)); 343 } 344 PetscCall(VecGetDM(X, &dmX)); 345 if (!dmX) dmX = dm; 346 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject*) §ion)); 347 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 348 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section,vStart,&bs)); 349 PetscCall(PetscSectionGetNumFields(section,&nfields)); 350 field = 0; 351 if (link->field >= 0) { 352 field = link->field; 353 nfields = field + 1; 354 } 355 for (i=0; field<(nfields?nfields:1); field++) { 356 PetscInt fbs,j; 357 const char *fieldname = NULL; 358 char buf[256]; 359 if (nfields) { /* We have user-defined fields/components */ 360 PetscCall(PetscSectionGetFieldDof(section,vStart,field,&fbs)); 361 PetscCall(PetscSectionGetFieldName(section,field,&fieldname)); 362 } else fbs = bs; /* Say we have one field with 'bs' components */ 363 if (nfields && !fieldname) { 364 PetscCall(PetscSNPrintf(buf,sizeof(buf),"PointField%" PetscInt_FMT,field)); 365 fieldname = buf; 366 } 367 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 368 PetscCheck(fbs <= 3,PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Point vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs); 369 #if defined(PETSC_USE_COMPLEX) 370 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,vecname,fieldname,boffset)); 371 boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int); 372 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,vecname,fieldname,boffset)); 373 boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int); 374 #else 375 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,vecname,fieldname,boffset)); 376 boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int); 377 #endif 378 } else { 379 for (j=0; j<fbs; j++) { 380 const char *compName = NULL; 381 char finalname[256]; 382 PetscCall(PetscSectionGetComponentName(section,field,j,&compName)); 383 PetscCall(PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName)); 384 #if defined(PETSC_USE_COMPLEX) 385 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,finalname,boffset)); 386 boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int); 387 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,finalname,boffset)); 388 boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int); 389 #else 390 PetscCall(PetscFPrintf(comm,fp," <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",precision,finalname,boffset)); 391 boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int); 392 #endif 393 } 394 } 395 } 396 } 397 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," </PointData>\n")); 398 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fp," </Piece>\n")); 399 } 400 } 401 402 PetscCall(PetscFPrintf(comm,fp," </UnstructuredGrid>\n")); 403 PetscCall(PetscFPrintf(comm,fp," <AppendedData encoding=\"raw\">\n")); 404 PetscCall(PetscFPrintf(comm,fp,"_")); 405 406 if (rank == 0) { 407 PetscInt maxsize = 0; 408 for (r=0; r<size; r++) { 409 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscVTUReal))); 410 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*3*sizeof(PetscVTUReal))); 411 maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt))); 412 } 413 PetscCall(PetscMalloc(maxsize,&buffer)); 414 } 415 for (r=0; r<size; r++) { 416 if (r == rank) { 417 PetscInt nsend; 418 { /* Position */ 419 const PetscScalar *x, *cx = NULL; 420 PetscVTUReal *y = NULL; 421 Vec coords, cellCoords; 422 PetscBool copy; 423 424 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 425 PetscCall(VecGetArrayRead(coords,&x)); 426 PetscCall(DMGetCellCoordinatesLocal(dm,&cellCoords)); 427 if (cellCoords) PetscCall(VecGetArrayRead(cellCoords,&cx)); 428 #if defined(PETSC_USE_COMPLEX) 429 copy = PETSC_TRUE; 430 #else 431 copy = (PetscBool) (dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal))); 432 #endif 433 if (copy) { 434 PetscCall(PetscMalloc1(piece.nvertices*3,&y)); 435 if (localized) { 436 PetscInt cnt; 437 for (c=cStart,cnt=0; c<cEnd; c++) { 438 PetscInt off, dof; 439 440 PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof)); 441 if (!dof) { 442 PetscInt *closure = NULL; 443 PetscInt closureSize; 444 445 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 446 for (v = 0; v < closureSize*2; v += 2) { 447 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 448 PetscCall(PetscSectionGetOffset(coordSection, closure[v], &off)); 449 if (dimEmbed != 3) { 450 y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]); 451 y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off+1]) : 0.0); 452 y[cnt*3+2] = (PetscVTUReal) 0.0; 453 } else { 454 y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]); 455 y[cnt*3+1] = (PetscVTUReal) PetscRealPart(x[off+1]); 456 y[cnt*3+2] = (PetscVTUReal) PetscRealPart(x[off+2]); 457 } 458 cnt++; 459 } 460 } 461 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure)); 462 } else { 463 PetscCall(PetscSectionGetOffset(cellCoordSection, c, &off)); 464 if (dimEmbed != 3) { 465 for (i=0; i<dof/dimEmbed; i++) { 466 y[cnt*3+0] = (PetscVTUReal) PetscRealPart(cx[off + i*dimEmbed + 0]); 467 y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(cx[off + i*dimEmbed + 1]) : 0.0); 468 y[cnt*3+2] = (PetscVTUReal) 0.0; 469 cnt++; 470 } 471 } else { 472 for (i=0; i<dof; i ++) { 473 y[cnt*3+i] = (PetscVTUReal) PetscRealPart(cx[off + i]); 474 } 475 cnt += dof/dimEmbed; 476 } 477 } 478 } 479 PetscCheck(cnt == piece.nvertices,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 480 } else { 481 for (i=0; i<piece.nvertices; i++) { 482 y[i*3+0] = (PetscVTUReal) PetscRealPart(x[i*dimEmbed+0]); 483 y[i*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[i*dimEmbed+1]) : 0.); 484 y[i*3+2] = (PetscVTUReal) ((dimEmbed > 2) ? PetscRealPart(x[i*dimEmbed+2]) : 0.); 485 } 486 } 487 } 488 nsend = piece.nvertices*3; 489 PetscCall(TransferWrite(comm,viewer,fp,r,0,copy ? (const void *) y : (const void *) x,buffer,nsend,MPIU_VTUREAL,tag)); 490 PetscCall(PetscFree(y)); 491 PetscCall(VecRestoreArrayRead(coords,&x)); 492 if (cellCoords) PetscCall(VecRestoreArrayRead(cellCoords,&cx)); 493 } 494 { /* Connectivity, offsets, types */ 495 PetscVTKInt *connectivity = NULL, *offsets = NULL; 496 PetscVTKType *types = NULL; 497 PetscCall(DMPlexGetVTKConnectivity(dm,localized,&piece,&connectivity,&offsets,&types)); 498 PetscCall(TransferWrite(comm,viewer,fp,r,0,connectivity,buffer,piece.nconn,MPI_INT,tag)); 499 PetscCall(TransferWrite(comm,viewer,fp,r,0,offsets,buffer,piece.ncells,MPI_INT,tag)); 500 PetscCall(TransferWrite(comm,viewer,fp,r,0,types,buffer,piece.ncells,MPI_CHAR,tag)); 501 PetscCall(PetscFree3(connectivity,offsets,types)); 502 } 503 { /* Owners (cell data) */ 504 PetscVTKInt *owners; 505 PetscCall(PetscMalloc1(piece.ncells,&owners)); 506 for (i=0; i<piece.ncells; i++) owners[i] = rank; 507 PetscCall(TransferWrite(comm,viewer,fp,r,0,owners,buffer,piece.ncells,MPI_INT,tag)); 508 PetscCall(PetscFree(owners)); 509 } 510 /* Cell data */ 511 for (link=vtk->link; link; link=link->next) { 512 Vec X = (Vec)link->vec; 513 DM dmX; 514 const PetscScalar *x; 515 PetscVTUReal *y; 516 PetscInt bs = 1, nfields, field; 517 PetscSection section = NULL; 518 519 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 520 PetscCall(VecGetDM(X, &dmX)); 521 if (!dmX) dmX = dm; 522 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject*) §ion)); 523 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 524 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section,cStart,&bs)); 525 PetscCall(PetscSectionGetNumFields(section,&nfields)); 526 field = 0; 527 if (link->field >= 0) { 528 field = link->field; 529 nfields = field + 1; 530 } 531 PetscCall(VecGetArrayRead(X,&x)); 532 PetscCall(PetscMalloc1(piece.ncells*3,&y)); 533 for (i=0; field<(nfields?nfields:1); field++) { 534 PetscInt fbs,j; 535 PetscFV fv = NULL; 536 PetscObject f; 537 PetscClassId fClass; 538 PetscBool vector; 539 if (nfields && cEnd > cStart) { /* We have user-defined fields/components */ 540 PetscCall(PetscSectionGetFieldDof(section,cStart,field,&fbs)); 541 } else fbs = bs; /* Say we have one field with 'bs' components */ 542 PetscCall(DMGetField(dmX,field,NULL,&f)); 543 PetscCall(PetscObjectGetClassId(f,&fClass)); 544 if (fClass == PETSCFV_CLASSID) { 545 fv = (PetscFV) f; 546 } 547 vector = PETSC_FALSE; 548 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 549 vector = PETSC_TRUE; 550 for (j = 0; j < fbs; j++) { 551 const char *compName = NULL; 552 if (fv) { 553 PetscCall(PetscFVGetComponentName(fv,j,&compName)); 554 if (compName) break; 555 } 556 } 557 if (j < fbs) vector = PETSC_FALSE; 558 } 559 if (vector) { 560 PetscInt cnt, l; 561 for (l = 0; l < loops_per_scalar; l++) { 562 for (c=cStart,cnt=0; c<cEnd; c++) { 563 const PetscScalar *xpoint; 564 PetscInt off, j; 565 566 if (hasLabel) { /* Ignore some cells */ 567 PetscInt value; 568 PetscCall(DMGetLabelValue(dmX, "vtk", c, &value)); 569 if (value != 1) continue; 570 } 571 if (nfields) { 572 PetscCall(PetscSectionGetFieldOffset(section, c, field, &off)); 573 } else { 574 PetscCall(PetscSectionGetOffset(section, c, &off)); 575 } 576 xpoint = &x[off]; 577 for (j = 0; j < fbs; j++) { 578 y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j])); 579 } 580 for (; j < 3; j++) y[cnt++] = 0.; 581 } 582 PetscCheck(cnt == piece.ncells*3,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 583 PetscCall(TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.ncells*3,MPIU_VTUREAL,tag)); 584 } 585 } else { 586 for (i=0; i<fbs; i++) { 587 PetscInt cnt, l; 588 for (l = 0; l < loops_per_scalar; l++) { 589 for (c=cStart,cnt=0; c<cEnd; c++) { 590 const PetscScalar *xpoint; 591 PetscInt off; 592 593 if (hasLabel) { /* Ignore some cells */ 594 PetscInt value; 595 PetscCall(DMGetLabelValue(dmX, "vtk", c, &value)); 596 if (value != 1) continue; 597 } 598 if (nfields) { 599 PetscCall(PetscSectionGetFieldOffset(section, c, field, &off)); 600 } else { 601 PetscCall(PetscSectionGetOffset(section, c, &off)); 602 } 603 xpoint = &x[off]; 604 y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i])); 605 } 606 PetscCheck(cnt == piece.ncells,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 607 PetscCall(TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.ncells,MPIU_VTUREAL,tag)); 608 } 609 } 610 } 611 } 612 PetscCall(PetscFree(y)); 613 PetscCall(VecRestoreArrayRead(X,&x)); 614 } 615 /* point data */ 616 for (link=vtk->link; link; link=link->next) { 617 Vec X = (Vec)link->vec; 618 DM dmX; 619 const PetscScalar *x; 620 PetscVTUReal *y; 621 PetscInt bs = 1, nfields, field; 622 PetscSection section = NULL; 623 624 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 625 PetscCall(VecGetDM(X, &dmX)); 626 if (!dmX) dmX = dm; 627 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject*) §ion)); 628 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 629 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section,vStart,&bs)); 630 PetscCall(PetscSectionGetNumFields(section,&nfields)); 631 field = 0; 632 if (link->field >= 0) { 633 field = link->field; 634 nfields = field + 1; 635 } 636 PetscCall(VecGetArrayRead(X,&x)); 637 PetscCall(PetscMalloc1(piece.nvertices*3,&y)); 638 for (i=0; field<(nfields?nfields:1); field++) { 639 PetscInt fbs,j; 640 if (nfields && vEnd > vStart) { /* We have user-defined fields/components */ 641 PetscCall(PetscSectionGetFieldDof(section,vStart,field,&fbs)); 642 } else fbs = bs; /* Say we have one field with 'bs' components */ 643 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 644 PetscInt cnt, l; 645 for (l = 0; l < loops_per_scalar; l++) { 646 if (!localized) { 647 for (v=vStart,cnt=0; v<vEnd; v++) { 648 PetscInt off; 649 const PetscScalar *xpoint; 650 651 if (nfields) { 652 PetscCall(PetscSectionGetFieldOffset(section,v,field,&off)); 653 } else { 654 PetscCall(PetscSectionGetOffset(section,v,&off)); 655 } 656 xpoint = &x[off]; 657 for (j = 0; j < fbs; j++) { 658 y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j])); 659 } 660 for (; j < 3; j++) y[cnt++] = 0.; 661 } 662 } else { 663 for (c=cStart,cnt=0; c<cEnd; c++) { 664 PetscInt *closure = NULL; 665 PetscInt closureSize, off; 666 667 PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 668 for (v = 0, off = 0; v < closureSize*2; v += 2) { 669 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 670 PetscInt voff; 671 const PetscScalar *xpoint; 672 673 if (nfields) { 674 PetscCall(PetscSectionGetFieldOffset(section,closure[v],field,&voff)); 675 } else { 676 PetscCall(PetscSectionGetOffset(section,closure[v],&voff)); 677 } 678 xpoint = &x[voff]; 679 for (j = 0; j < fbs; j++) { 680 y[cnt + off++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j])); 681 } 682 for (; j < 3; j++) y[cnt + off++] = 0.; 683 } 684 } 685 cnt += off; 686 PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 687 } 688 } 689 PetscCheck(cnt == piece.nvertices*3,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 690 PetscCall(TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.nvertices*3,MPIU_VTUREAL,tag)); 691 } 692 } else { 693 for (i=0; i<fbs; i++) { 694 PetscInt cnt, l; 695 for (l = 0; l < loops_per_scalar; l++) { 696 if (!localized) { 697 for (v=vStart,cnt=0; v<vEnd; v++) { 698 PetscInt off; 699 const PetscScalar *xpoint; 700 701 if (nfields) { 702 PetscCall(PetscSectionGetFieldOffset(section,v,field,&off)); 703 } else { 704 PetscCall(PetscSectionGetOffset(section,v,&off)); 705 } 706 xpoint = &x[off]; 707 y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i])); 708 } 709 } else { 710 for (c=cStart,cnt=0; c<cEnd; c++) { 711 PetscInt *closure = NULL; 712 PetscInt closureSize, off; 713 714 PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 715 for (v = 0, off = 0; v < closureSize*2; v += 2) { 716 if ((closure[v] >= vStart) && (closure[v] < vEnd)) { 717 PetscInt voff; 718 const PetscScalar *xpoint; 719 720 if (nfields) { 721 PetscCall(PetscSectionGetFieldOffset(section,closure[v],field,&voff)); 722 } else { 723 PetscCall(PetscSectionGetOffset(section,closure[v],&voff)); 724 } 725 xpoint = &x[voff]; 726 y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i])); 727 } 728 } 729 cnt += off; 730 PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure)); 731 } 732 } 733 PetscCheck(cnt == piece.nvertices,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match"); 734 PetscCall(TransferWrite(comm,viewer,fp,r,0,y,buffer,piece.nvertices,MPIU_VTUREAL,tag)); 735 } 736 } 737 } 738 } 739 PetscCall(PetscFree(y)); 740 PetscCall(VecRestoreArrayRead(X,&x)); 741 } 742 } else if (rank == 0) { 743 PetscInt l; 744 745 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag)); /* positions */ 746 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,MPI_INT,tag)); /* connectivity */ 747 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag)); /* offsets */ 748 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_CHAR,tag)); /* types */ 749 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag)); /* owner rank (cells) */ 750 /* all cell data */ 751 for (link=vtk->link; link; link=link->next) { 752 Vec X = (Vec)link->vec; 753 PetscInt bs = 1, nfields, field; 754 DM dmX; 755 PetscSection section = NULL; 756 757 if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue; 758 PetscCall(VecGetDM(X, &dmX)); 759 if (!dmX) dmX = dm; 760 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject*) §ion)); 761 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 762 if (cEnd > cStart) PetscCall(PetscSectionGetDof(section,cStart,&bs)); 763 PetscCall(PetscSectionGetNumFields(section,&nfields)); 764 field = 0; 765 if (link->field >= 0) { 766 field = link->field; 767 nfields = field + 1; 768 } 769 for (i=0; field<(nfields?nfields:1); field++) { 770 PetscInt fbs,j; 771 PetscFV fv = NULL; 772 PetscObject f; 773 PetscClassId fClass; 774 PetscBool vector; 775 if (nfields && cEnd > cStart) { /* We have user-defined fields/components */ 776 PetscCall(PetscSectionGetFieldDof(section,cStart,field,&fbs)); 777 } else fbs = bs; /* Say we have one field with 'bs' components */ 778 PetscCall(DMGetField(dmX,field,NULL,&f)); 779 PetscCall(PetscObjectGetClassId(f,&fClass)); 780 if (fClass == PETSCFV_CLASSID) { 781 fv = (PetscFV) f; 782 } 783 vector = PETSC_FALSE; 784 if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) { 785 vector = PETSC_TRUE; 786 for (j = 0; j < fbs; j++) { 787 const char *compName = NULL; 788 if (fv) { 789 PetscCall(PetscFVGetComponentName(fv,j,&compName)); 790 if (compName) break; 791 } 792 } 793 if (j < fbs) vector = PETSC_FALSE; 794 } 795 if (vector) { 796 for (l = 0; l < loops_per_scalar; l++) { 797 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells*3,MPIU_VTUREAL,tag)); 798 } 799 } else { 800 for (i=0; i<fbs; i++) { 801 for (l = 0; l < loops_per_scalar; l++) { 802 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPIU_VTUREAL,tag)); 803 } 804 } 805 } 806 } 807 } 808 /* all point data */ 809 for (link=vtk->link; link; link=link->next) { 810 Vec X = (Vec)link->vec; 811 DM dmX; 812 PetscInt bs = 1, nfields, field; 813 PetscSection section = NULL; 814 815 if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue; 816 PetscCall(VecGetDM(X, &dmX)); 817 if (!dmX) dmX = dm; 818 PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject*) §ion)); 819 if (!section) PetscCall(DMGetLocalSection(dmX, §ion)); 820 if (vEnd > vStart) PetscCall(PetscSectionGetDof(section,vStart,&bs)); 821 PetscCall(PetscSectionGetNumFields(section,&nfields)); 822 field = 0; 823 if (link->field >= 0) { 824 field = link->field; 825 nfields = field + 1; 826 } 827 for (i=0; field<(nfields?nfields:1); field++) { 828 PetscInt fbs; 829 if (nfields && vEnd > vStart) { /* We have user-defined fields/components */ 830 PetscCall(PetscSectionGetFieldDof(section,vStart,field,&fbs)); 831 } else fbs = bs; /* Say we have one field with 'bs' components */ 832 if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) { 833 for (l = 0; l < loops_per_scalar; l++) { 834 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag)); 835 } 836 } else { 837 for (i=0; i<fbs; i++) { 838 for (l = 0; l < loops_per_scalar; l++) { 839 PetscCall(TransferWrite(comm,viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,MPIU_VTUREAL,tag)); 840 } 841 } 842 } 843 } 844 } 845 } 846 } 847 PetscCall(PetscFree(gpiece)); 848 PetscCall(PetscFree(buffer)); 849 PetscCall(PetscFPrintf(comm,fp,"\n </AppendedData>\n")); 850 PetscCall(PetscFPrintf(comm,fp,"</VTKFile>\n")); 851 PetscCall(PetscFClose(comm,fp)); 852 PetscFunctionReturn(0); 853 } 854