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