Actual source code: plexvtu.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

      4: typedef struct {
      5:   PetscInt nvertices;
      6:   PetscInt ncells;
      7:   PetscInt nconn; /* number of entries in cell->vertex connectivity array */
      8: } PieceInfo;

     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

     26: static PetscErrorCode TransferWrite(MPI_Comm comm, PetscViewer viewer, FILE *fp, PetscMPIInt srank, PetscMPIInt root, const void *send, void *recv, PetscCount count, MPI_Datatype mpidatatype, PetscMPIInt tag)
     27: {
     28:   PetscMPIInt rank;

     30:   PetscFunctionBegin;
     31:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     32:   if (rank == srank && rank != root) {
     33:     PetscCallMPI(MPIU_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(MPIU_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(PETSC_SUCCESS);
     49: }

     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;

     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;

     70:   countcell = 0;
     71:   countconn = 0;
     72:   for (c = cStart; c < cEnd; ++c) {
     73:     PetscInt nverts, dof = 0, celltype, startoffset, nC = 0;

     75:     if (hasLabel) {
     76:       PetscInt value;

     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;

     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) PetscCall(PetscVTKIntCast(closure[v] - vStart, &conn[countconn++]));
     91:           else PetscCall(PetscVTKIntCast(startoffset + nC, &conn[countconn++]));
     92:           ++nC;
     93:         }
     94:       }
     95:       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
     96:     } else {
     97:       for (nC = 0; nC < dof / dim; nC++) PetscCall(PetscVTKIntCast(startoffset + nC, &conn[countconn++]));
     98:     }

    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) PetscCall(PetscVTKIntCast(cone[i], &conn[s + i]));
    105:     }
    106:     PetscCall(PetscVTKIntCast(countconn, &offsets[countcell]));

    108:     nverts = countconn - startoffset;
    109:     PetscCall(DMPlexVTKGetCellType_Internal(dm, dim, nverts, &celltype));

    111:     types[countcell] = (PetscVTKType)celltype;
    112:     countcell++;
    113:   }
    114:   PetscCheck(countcell == piece->ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent cell count");
    115:   PetscCheck(countconn == piece->nconn, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent connectivity count");
    116:   *oconn    = conn;
    117:   *ooffsets = offsets;
    118:   *otypes   = types;
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: PETSC_INTERN PetscErrorCode DMPlexGetNonEmptyComm_Private(DM dm, MPI_Comm *comm)
    123: {
    124:   DM_Plex *mesh = (DM_Plex *)dm->data;

    126:   PetscFunctionBegin;
    127:   if (mesh->nonempty_comm == MPI_COMM_SELF) { /* Not yet setup */
    128:     PetscInt    cStart, cEnd, cellHeight;
    129:     MPI_Comm    dmcomm = PetscObjectComm((PetscObject)dm);
    130:     PetscMPIInt color, rank;

    132:     PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
    133:     PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
    134:     color = (cStart < cEnd) ? 0 : 1;
    135:     PetscCallMPI(MPI_Comm_rank(dmcomm, &rank));
    136:     PetscCallMPI(MPI_Comm_split(dmcomm, color, rank, &mesh->nonempty_comm));
    137:     if (color == 1) {
    138:       PetscCallMPI(MPI_Comm_free(&mesh->nonempty_comm));
    139:       mesh->nonempty_comm = MPI_COMM_NULL;
    140:     }
    141:   }
    142:   *comm = mesh->nonempty_comm;
    143:   PetscFunctionReturn(PETSC_SUCCESS);
    144: }

    146: static PetscErrorCode DMGetFieldIfFV_Private(DM dm, PetscInt field, PetscFV *fv)
    147: {
    148:   PetscObject  f      = NULL;
    149:   PetscClassId fClass = PETSC_SMALLEST_CLASSID;
    150:   PetscInt     nf;

    152:   PetscFunctionBegin;
    153:   *fv = NULL;
    154:   PetscCall(DMGetNumFields(dm, &nf));
    155:   if (nf > 0) {
    156:     PetscCall(DMGetField(dm, field, NULL, &f));
    157:     PetscCall(PetscObjectGetClassId(f, &fClass));
    158:     if (fClass == PETSCFV_CLASSID) *fv = (PetscFV)f;
    159:   }
    160:   PetscFunctionReturn(PETSC_SUCCESS);
    161: }

    163: /*
    164:   Write all fields that have been provided to the viewer
    165:   Multi-block XML format with binary appended data.
    166: */
    167: PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm, PetscViewer viewer)
    168: {
    169:   MPI_Comm                 comm;
    170:   PetscSection             coordSection, cellCoordSection;
    171:   PetscViewer_VTK         *vtk = (PetscViewer_VTK *)viewer->data;
    172:   PetscViewerVTKObjectLink link;
    173:   FILE                    *fp;
    174:   PetscMPIInt              rank, size, tag;
    175:   PetscInt                 dimEmbed, cellHeight, cStart, cEnd, vStart, vEnd, numLabelCells, hasLabel, c, v, i;
    176:   PetscBool                localized;
    177:   PieceInfo                piece, *gpiece = NULL;
    178:   void                    *buffer     = NULL;
    179:   const char              *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
    180:   PetscInt                 loops_per_scalar;

    182:   PetscFunctionBegin;
    183:   PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
    184:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
    185:   PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
    186:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    187:   PetscCall(DMGetStratumSize(dm, "vtk", 1, &numLabelCells));
    188:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
    189:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
    190:   PetscCall(DMGetCellCoordinateSection(dm, &cellCoordSection));
    191:   PetscCall(PetscCommGetNewTag(PetscObjectComm((PetscObject)dm), &tag));

    193:   PetscCall(DMPlexGetNonEmptyComm_Private(dm, &comm));
    194: #if defined(PETSC_USE_COMPLEX)
    195:   loops_per_scalar = 2;
    196: #else
    197:   loops_per_scalar = 1;
    198: #endif
    199:   if (comm == MPI_COMM_NULL) goto finalize;
    200:   PetscCallMPI(MPI_Comm_size(comm, &size));
    201:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    203:   PetscCall(PetscFOpen(comm, vtk->filename, "wb", &fp));
    204:   PetscCall(PetscFPrintf(comm, fp, "<?xml version=\"1.0\"?>\n"));
    205:   PetscCall(PetscFPrintf(comm, fp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\" header_type=\"UInt64\">\n", byte_order));
    206:   PetscCall(PetscFPrintf(comm, fp, "  <UnstructuredGrid>\n"));

    208:   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
    209:   piece.nvertices = 0;
    210:   piece.ncells    = 0;
    211:   piece.nconn     = 0;
    212:   if (!localized) piece.nvertices = vEnd - vStart;
    213:   for (c = cStart; c < cEnd; ++c) {
    214:     PetscInt dof = 0;
    215:     if (hasLabel) {
    216:       PetscInt value;

    218:       PetscCall(DMGetLabelValue(dm, "vtk", c, &value));
    219:       if (value != 1) continue;
    220:     }
    221:     if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
    222:     if (!dof) {
    223:       PetscInt *closure = NULL;
    224:       PetscInt  closureSize;

    226:       PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    227:       for (v = 0; v < closureSize * 2; v += 2) {
    228:         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
    229:           piece.nconn++;
    230:           if (localized) piece.nvertices++;
    231:         }
    232:       }
    233:       PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    234:     } else {
    235:       piece.nvertices += dof / dimEmbed;
    236:       piece.nconn += dof / dimEmbed;
    237:     }
    238:     piece.ncells++;
    239:   }
    240:   if (rank == 0) PetscCall(PetscMalloc1(size, &gpiece));
    241:   PetscCallMPI(MPI_Gather((PetscInt *)&piece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, (PetscInt *)gpiece, sizeof(piece) / sizeof(PetscInt), MPIU_INT, 0, comm));

    243:   /*
    244:    * Write file header
    245:    */
    246:   if (rank == 0) {
    247:     PetscInt64 boffset = 0;

    249:     for (PetscMPIInt r = 0; r < size; r++) {
    250:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "    <Piece NumberOfPoints=\"%" PetscInt_FMT "\" NumberOfCells=\"%" PetscInt_FMT "\">\n", gpiece[r].nvertices, gpiece[r].ncells));
    251:       /* Coordinate positions */
    252:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <Points>\n"));
    253:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, boffset));
    254:       boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    255:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </Points>\n"));
    256:       /* Cell connectivity */
    257:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <Cells>\n"));
    258:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
    259:       boffset += gpiece[r].nconn * sizeof(PetscVTKInt) + (gpiece[r].nconn ? sizeof(PetscInt64) : 0);
    260:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"offsets\"      NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
    261:       boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    262:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"UInt8\" Name=\"types\"        NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
    263:       boffset += gpiece[r].ncells * sizeof(unsigned char) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    264:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </Cells>\n"));

    266:       /*
    267:        * Cell Data headers
    268:        */
    269:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <CellData>\n"));
    270:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", boffset));
    271:       boffset += gpiece[r].ncells * sizeof(PetscVTKInt) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    272:       /* all the vectors */
    273:       for (link = vtk->link; link; link = link->next) {
    274:         Vec          X       = (Vec)link->vec;
    275:         DM           dmX     = NULL;
    276:         PetscInt     bs      = 1, nfields, field;
    277:         const char  *vecname = "";
    278:         PetscSection section;
    279:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
    280:         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. */
    281:           PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
    282:         }
    283:         PetscCall(VecGetDM(X, &dmX));
    284:         if (!dmX) dmX = dm;
    285:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    286:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    287:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
    288:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    289:         field = 0;
    290:         if (link->field >= 0) {
    291:           field   = link->field;
    292:           nfields = field + 1;
    293:         }
    294:         for (i = 0; field < (nfields ? nfields : 1); field++) {
    295:           PetscInt    fbs, j;
    296:           PetscFV     fv        = NULL;
    297:           const char *fieldname = NULL;
    298:           char        buf[256];
    299:           PetscBool   vector;

    301:           if (nfields) { /* We have user-defined fields/components */
    302:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
    303:             PetscCall(PetscSectionGetFieldName(section, field, &fieldname));
    304:           } else fbs = bs; /* Say we have one field with 'bs' components */
    305:           PetscCall(DMGetFieldIfFV_Private(dmX, field, &fv));
    306:           if (nfields && !fieldname) {
    307:             PetscCall(PetscSNPrintf(buf, sizeof(buf), "CellField%" PetscInt_FMT, field));
    308:             fieldname = buf;
    309:           }
    310:           vector = PETSC_FALSE;
    311:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
    312:             vector = PETSC_TRUE;
    313:             PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Cell vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs);
    314:             for (j = 0; j < fbs; j++) {
    315:               const char *compName = NULL;
    316:               if (fv) {
    317:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
    318:                 if (compName) break;
    319:               }
    320:             }
    321:             if (j < fbs) vector = PETSC_FALSE;
    322:           }
    323:           if (vector) {
    324: #if defined(PETSC_USE_COMPLEX)
    325:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    326:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    327:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    328:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    329: #else
    330:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    331:             boffset += gpiece[r].ncells * 3 * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    332: #endif
    333:             i += fbs;
    334:           } else {
    335:             for (j = 0; j < fbs; j++) {
    336:               const char *compName = NULL;
    337:               char        finalname[256];
    338:               if (fv) PetscCall(PetscFVGetComponentName(fv, j, &compName));
    339:               if (compName) {
    340:                 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName));
    341:               } else if (fbs > 1) {
    342:                 PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%" PetscInt_FMT, vecname, fieldname, j));
    343:               } else {
    344:                 PetscCall(PetscSNPrintf(finalname, 255, "%s%s", vecname, fieldname));
    345:               }
    346: #if defined(PETSC_USE_COMPLEX)
    347:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
    348:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    349:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
    350:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    351: #else
    352:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
    353:               boffset += gpiece[r].ncells * sizeof(PetscVTUReal) + (gpiece[r].ncells ? sizeof(PetscInt64) : 0);
    354: #endif
    355:               i++;
    356:             }
    357:           }
    358:         }
    359:         //PetscCheck(i == bs,comm,PETSC_ERR_PLIB,"Total number of field components %" PetscInt_FMT " != block size %" PetscInt_FMT,i,bs);
    360:       }
    361:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </CellData>\n"));

    363:       /*
    364:        * Point Data headers
    365:        */
    366:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      <PointData>\n"));
    367:       for (link = vtk->link; link; link = link->next) {
    368:         Vec          X = (Vec)link->vec;
    369:         DM           dmX;
    370:         PetscInt     bs      = 1, nfields, field;
    371:         const char  *vecname = "";
    372:         PetscSection section;
    373:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
    374:         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. */
    375:           PetscCall(PetscObjectGetName((PetscObject)X, &vecname));
    376:         }
    377:         PetscCall(VecGetDM(X, &dmX));
    378:         if (!dmX) dmX = dm;
    379:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    380:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    381:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
    382:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    383:         field = 0;
    384:         if (link->field >= 0) {
    385:           field   = link->field;
    386:           nfields = field + 1;
    387:         }
    388:         for (; field < (nfields ? nfields : 1); field++) {
    389:           PetscInt    fbs, j;
    390:           const char *fieldname = NULL;
    391:           char        buf[256];
    392:           if (nfields) { /* We have user-defined fields/components */
    393:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
    394:             PetscCall(PetscSectionGetFieldName(section, field, &fieldname));
    395:           } else fbs = bs; /* Say we have one field with 'bs' components */
    396:           if (nfields && !fieldname) {
    397:             PetscCall(PetscSNPrintf(buf, sizeof(buf), "PointField%" PetscInt_FMT, field));
    398:             fieldname = buf;
    399:           }
    400:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
    401:             PetscCheck(fbs <= 3, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_SIZ, "Point vector fields can have at most 3 components, %" PetscInt_FMT " given", fbs);
    402: #if defined(PETSC_USE_COMPLEX)
    403:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    404:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    405:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    406:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    407: #else
    408:             PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, vecname, fieldname, boffset));
    409:             boffset += gpiece[r].nvertices * 3 * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    410: #endif
    411:           } else {
    412:             for (j = 0; j < fbs; j++) {
    413:               const char *compName = NULL;
    414:               char        finalname[256];
    415:               PetscCall(PetscSectionGetComponentName(section, field, j, &compName));
    416:               PetscCall(PetscSNPrintf(finalname, 255, "%s%s.%s", vecname, fieldname, compName));
    417: #if defined(PETSC_USE_COMPLEX)
    418:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
    419:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    420:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
    421:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    422: #else
    423:               PetscCall(PetscFPrintf(comm, fp, "        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%" PetscInt64_FMT "\" />\n", precision, finalname, boffset));
    424:               boffset += gpiece[r].nvertices * sizeof(PetscVTUReal) + (gpiece[r].nvertices ? sizeof(PetscInt64) : 0);
    425: #endif
    426:             }
    427:           }
    428:         }
    429:       }
    430:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "      </PointData>\n"));
    431:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "    </Piece>\n"));
    432:     }
    433:   }

    435:   PetscCall(PetscFPrintf(comm, fp, "  </UnstructuredGrid>\n"));
    436:   PetscCall(PetscFPrintf(comm, fp, "  <AppendedData encoding=\"raw\">\n"));
    437:   PetscCall(PetscFPrintf(comm, fp, "_"));

    439:   if (rank == 0) {
    440:     PetscInt maxsize = 0;
    441:     for (PetscMPIInt r = 0; r < size; r++) {
    442:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal)));
    443:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal)));
    444:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt)));
    445:     }
    446:     PetscCall(PetscMalloc(maxsize, &buffer));
    447:   }
    448:   for (PetscMPIInt r = 0; r < size; r++) {
    449:     if (r == rank) {
    450:       PetscInt nsend;
    451:       { /* Position */
    452:         const PetscScalar *x, *cx = NULL;
    453:         PetscVTUReal      *y = NULL;
    454:         Vec                coords, cellCoords;
    455:         PetscBool          copy;

    457:         PetscCall(DMGetCoordinatesLocal(dm, &coords));
    458:         PetscCall(VecGetArrayRead(coords, &x));
    459:         PetscCall(DMGetCellCoordinatesLocal(dm, &cellCoords));
    460:         if (cellCoords) PetscCall(VecGetArrayRead(cellCoords, &cx));
    461: #if defined(PETSC_USE_COMPLEX)
    462:         copy = PETSC_TRUE;
    463: #else
    464:         copy = (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
    465: #endif
    466:         if (copy) {
    467:           PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
    468:           if (localized) {
    469:             PetscInt cnt;
    470:             for (c = cStart, cnt = 0; c < cEnd; c++) {
    471:               PetscInt off, dof;

    473:               PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
    474:               if (!dof) {
    475:                 PetscInt *closure = NULL;
    476:                 PetscInt  closureSize;

    478:                 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    479:                 for (v = 0; v < closureSize * 2; v += 2) {
    480:                   if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
    481:                     PetscCall(PetscSectionGetOffset(coordSection, closure[v], &off));
    482:                     if (dimEmbed != 3) {
    483:                       y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
    484:                       y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[off + 1]) : 0.0);
    485:                       y[cnt * 3 + 2] = (PetscVTUReal)0.0;
    486:                     } else {
    487:                       y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(x[off + 0]);
    488:                       y[cnt * 3 + 1] = (PetscVTUReal)PetscRealPart(x[off + 1]);
    489:                       y[cnt * 3 + 2] = (PetscVTUReal)PetscRealPart(x[off + 2]);
    490:                     }
    491:                     cnt++;
    492:                   }
    493:                 }
    494:                 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
    495:               } else {
    496:                 PetscCall(PetscSectionGetOffset(cellCoordSection, c, &off));
    497:                 if (dimEmbed != 3) {
    498:                   for (i = 0; i < dof / dimEmbed; i++) {
    499:                     y[cnt * 3 + 0] = (PetscVTUReal)PetscRealPart(cx[off + i * dimEmbed + 0]);
    500:                     y[cnt * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(cx[off + i * dimEmbed + 1]) : 0.0);
    501:                     y[cnt * 3 + 2] = (PetscVTUReal)0.0;
    502:                     cnt++;
    503:                   }
    504:                 } else {
    505:                   for (i = 0; i < dof; i++) y[cnt * 3 + i] = (PetscVTUReal)PetscRealPart(cx[off + i]);
    506:                   cnt += dof / dimEmbed;
    507:                 }
    508:               }
    509:             }
    510:             PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    511:           } else {
    512:             for (i = 0; i < piece.nvertices; i++) {
    513:               y[i * 3 + 0] = (PetscVTUReal)PetscRealPart(x[i * dimEmbed + 0]);
    514:               y[i * 3 + 1] = (PetscVTUReal)((dimEmbed > 1) ? PetscRealPart(x[i * dimEmbed + 1]) : 0.);
    515:               y[i * 3 + 2] = (PetscVTUReal)((dimEmbed > 2) ? PetscRealPart(x[i * dimEmbed + 2]) : 0.);
    516:             }
    517:           }
    518:         }
    519:         nsend = piece.nvertices * 3;
    520:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, copy ? (const void *)y : (const void *)x, buffer, nsend, MPIU_VTUREAL, tag));
    521:         PetscCall(PetscFree(y));
    522:         PetscCall(VecRestoreArrayRead(coords, &x));
    523:         if (cellCoords) PetscCall(VecRestoreArrayRead(cellCoords, &cx));
    524:       }
    525:       { /* Connectivity, offsets, types */
    526:         PetscVTKInt  *connectivity = NULL, *offsets = NULL;
    527:         PetscVTKType *types = NULL;
    528:         PetscCall(DMPlexGetVTKConnectivity(dm, localized, &piece, &connectivity, &offsets, &types));
    529:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, connectivity, buffer, piece.nconn, MPI_INT, tag));
    530:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, offsets, buffer, piece.ncells, MPI_INT, tag));
    531:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, types, buffer, piece.ncells, MPI_CHAR, tag));
    532:         PetscCall(PetscFree3(connectivity, offsets, types));
    533:       }
    534:       { /* Owners (cell data) */
    535:         PetscVTKInt *owners;
    536:         PetscMPIInt  orank;

    538:         PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &orank));
    539:         PetscCall(PetscMalloc1(piece.ncells, &owners));
    540:         for (i = 0; i < piece.ncells; i++) owners[i] = orank;
    541:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag));
    542:         PetscCall(PetscFree(owners));
    543:       }
    544:       /* Cell data */
    545:       for (link = vtk->link; link; link = link->next) {
    546:         Vec                X = (Vec)link->vec;
    547:         DM                 dmX;
    548:         const PetscScalar *x;
    549:         PetscVTUReal      *y;
    550:         PetscInt           bs      = 1, nfields, field;
    551:         PetscSection       section = NULL;

    553:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
    554:         PetscCall(VecGetDM(X, &dmX));
    555:         if (!dmX) dmX = dm;
    556:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    557:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    558:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
    559:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    560:         field = 0;
    561:         if (link->field >= 0) {
    562:           field   = link->field;
    563:           nfields = field + 1;
    564:         }
    565:         PetscCall(VecGetArrayRead(X, &x));
    566:         PetscCall(PetscMalloc1(piece.ncells * 3, &y));
    567:         for (; field < (nfields ? nfields : 1); field++) {
    568:           PetscInt  fbs, j;
    569:           PetscFV   fv = NULL;
    570:           PetscBool vector;

    572:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
    573:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
    574:           } else fbs = bs; /* Say we have one field with 'bs' components */
    575:           PetscCall(DMGetFieldIfFV_Private(dmX, field, &fv));
    576:           vector = PETSC_FALSE;
    577:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
    578:             vector = PETSC_TRUE;
    579:             for (j = 0; j < fbs; j++) {
    580:               const char *compName = NULL;
    581:               if (fv) {
    582:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
    583:                 if (compName) break;
    584:               }
    585:             }
    586:             if (j < fbs) vector = PETSC_FALSE;
    587:           }
    588:           if (vector) {
    589:             PetscInt cnt, l;
    590:             for (l = 0; l < loops_per_scalar; l++) {
    591:               for (c = cStart, cnt = 0; c < cEnd; c++) {
    592:                 const PetscScalar *xpoint;
    593:                 PetscInt           off, j;

    595:                 if (hasLabel) { /* Ignore some cells */
    596:                   PetscInt value;
    597:                   PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
    598:                   if (value != 1) continue;
    599:                 }
    600:                 if (nfields) {
    601:                   PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
    602:                 } else {
    603:                   PetscCall(PetscSectionGetOffset(section, c, &off));
    604:                 }
    605:                 xpoint = &x[off];
    606:                 for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
    607:                 for (; j < 3; j++) y[cnt++] = 0.;
    608:               }
    609:               PetscCheck(cnt == piece.ncells * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    610:               PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag));
    611:             }
    612:           } else {
    613:             for (i = 0; i < fbs; i++) {
    614:               PetscInt cnt, l;
    615:               for (l = 0; l < loops_per_scalar; l++) {
    616:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
    617:                   const PetscScalar *xpoint;
    618:                   PetscInt           off;

    620:                   if (hasLabel) { /* Ignore some cells */
    621:                     PetscInt value;
    622:                     PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
    623:                     if (value != 1) continue;
    624:                   }
    625:                   if (nfields) {
    626:                     PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
    627:                   } else {
    628:                     PetscCall(PetscSectionGetOffset(section, c, &off));
    629:                   }
    630:                   xpoint   = &x[off];
    631:                   y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
    632:                 }
    633:                 PetscCheck(cnt == piece.ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    634:                 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag));
    635:               }
    636:             }
    637:           }
    638:         }
    639:         PetscCall(PetscFree(y));
    640:         PetscCall(VecRestoreArrayRead(X, &x));
    641:       }
    642:       /* point data */
    643:       for (link = vtk->link; link; link = link->next) {
    644:         Vec                X = (Vec)link->vec;
    645:         DM                 dmX;
    646:         const PetscScalar *x;
    647:         PetscVTUReal      *y;
    648:         PetscInt           bs      = 1, nfields, field;
    649:         PetscSection       section = NULL;

    651:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
    652:         PetscCall(VecGetDM(X, &dmX));
    653:         if (!dmX) dmX = dm;
    654:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    655:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    656:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
    657:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    658:         field = 0;
    659:         if (link->field >= 0) {
    660:           field   = link->field;
    661:           nfields = field + 1;
    662:         }
    663:         PetscCall(VecGetArrayRead(X, &x));
    664:         PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
    665:         for (; field < (nfields ? nfields : 1); field++) {
    666:           PetscInt fbs, j;
    667:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
    668:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
    669:           } else fbs = bs; /* Say we have one field with 'bs' components */
    670:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
    671:             PetscInt cnt, l;
    672:             for (l = 0; l < loops_per_scalar; l++) {
    673:               if (!localized) {
    674:                 for (v = vStart, cnt = 0; v < vEnd; v++) {
    675:                   PetscInt           off;
    676:                   const PetscScalar *xpoint;

    678:                   if (nfields) {
    679:                     PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
    680:                   } else {
    681:                     PetscCall(PetscSectionGetOffset(section, v, &off));
    682:                   }
    683:                   xpoint = &x[off];
    684:                   for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
    685:                   for (; j < 3; j++) y[cnt++] = 0.;
    686:                 }
    687:               } else {
    688:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
    689:                   PetscInt *closure = NULL;
    690:                   PetscInt  closureSize, off;

    692:                   PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    693:                   for (v = 0, off = 0; v < closureSize * 2; v += 2) {
    694:                     if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
    695:                       PetscInt           voff;
    696:                       const PetscScalar *xpoint;

    698:                       if (nfields) {
    699:                         PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
    700:                       } else {
    701:                         PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
    702:                       }
    703:                       xpoint = &x[voff];
    704:                       for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
    705:                       for (; j < 3; j++) y[cnt + off++] = 0.;
    706:                     }
    707:                   }
    708:                   cnt += off;
    709:                   PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    710:                 }
    711:               }
    712:               PetscCheck(cnt == piece.nvertices * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    713:               PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag));
    714:             }
    715:           } else {
    716:             for (i = 0; i < fbs; i++) {
    717:               PetscInt cnt, l;
    718:               for (l = 0; l < loops_per_scalar; l++) {
    719:                 if (!localized) {
    720:                   for (v = vStart, cnt = 0; v < vEnd; v++) {
    721:                     PetscInt           off;
    722:                     const PetscScalar *xpoint;

    724:                     if (nfields) {
    725:                       PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
    726:                     } else {
    727:                       PetscCall(PetscSectionGetOffset(section, v, &off));
    728:                     }
    729:                     xpoint   = &x[off];
    730:                     y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
    731:                   }
    732:                 } else {
    733:                   for (c = cStart, cnt = 0; c < cEnd; c++) {
    734:                     PetscInt *closure = NULL;
    735:                     PetscInt  closureSize, off;

    737:                     PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    738:                     for (v = 0, off = 0; v < closureSize * 2; v += 2) {
    739:                       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
    740:                         PetscInt           voff;
    741:                         const PetscScalar *xpoint;

    743:                         if (nfields) {
    744:                           PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
    745:                         } else {
    746:                           PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
    747:                         }
    748:                         xpoint         = &x[voff];
    749:                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
    750:                       }
    751:                     }
    752:                     cnt += off;
    753:                     PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    754:                   }
    755:                 }
    756:                 PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    757:                 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag));
    758:               }
    759:             }
    760:           }
    761:         }
    762:         PetscCall(PetscFree(y));
    763:         PetscCall(VecRestoreArrayRead(X, &x));
    764:       }
    765:     } else if (rank == 0) {
    766:       PetscInt l;

    768:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); /* positions */
    769:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag));              /* connectivity */
    770:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag));             /* offsets */
    771:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag));            /* types */
    772:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag));             /* owner rank (cells) */
    773:       /* all cell data */
    774:       for (link = vtk->link; link; link = link->next) {
    775:         Vec          X  = (Vec)link->vec;
    776:         PetscInt     bs = 1, nfields, field;
    777:         DM           dmX;
    778:         PetscSection section = NULL;

    780:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
    781:         PetscCall(VecGetDM(X, &dmX));
    782:         if (!dmX) dmX = dm;
    783:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    784:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    785:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
    786:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    787:         field = 0;
    788:         if (link->field >= 0) {
    789:           field   = link->field;
    790:           nfields = field + 1;
    791:         }
    792:         for (; field < (nfields ? nfields : 1); field++) {
    793:           PetscInt  fbs, j;
    794:           PetscFV   fv = NULL;
    795:           PetscBool vector;

    797:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
    798:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
    799:           } else fbs = bs; /* Say we have one field with 'bs' components */
    800:           PetscCall(DMGetFieldIfFV_Private(dmX, field, &fv));
    801:           vector = PETSC_FALSE;
    802:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
    803:             vector = PETSC_TRUE;
    804:             for (j = 0; j < fbs; j++) {
    805:               const char *compName = NULL;
    806:               if (fv) {
    807:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
    808:                 if (compName) break;
    809:               }
    810:             }
    811:             if (j < fbs) vector = PETSC_FALSE;
    812:           }
    813:           if (vector) {
    814:             for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag));
    815:           } else {
    816:             for (i = 0; i < fbs; i++) {
    817:               for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag));
    818:             }
    819:           }
    820:         }
    821:       }
    822:       /* all point data */
    823:       for (link = vtk->link; link; link = link->next) {
    824:         Vec          X = (Vec)link->vec;
    825:         DM           dmX;
    826:         PetscInt     bs      = 1, nfields, field;
    827:         PetscSection section = NULL;

    829:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
    830:         PetscCall(VecGetDM(X, &dmX));
    831:         if (!dmX) dmX = dm;
    832:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    833:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    834:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
    835:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    836:         field = 0;
    837:         if (link->field >= 0) {
    838:           field   = link->field;
    839:           nfields = field + 1;
    840:         }
    841:         for (; field < (nfields ? nfields : 1); field++) {
    842:           PetscInt fbs;
    843:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
    844:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
    845:           } else fbs = bs; /* Say we have one field with 'bs' components */
    846:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
    847:             for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag));
    848:           } else {
    849:             for (i = 0; i < fbs; i++) {
    850:               for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag));
    851:             }
    852:           }
    853:         }
    854:       }
    855:     }
    856:   }
    857:   PetscCall(PetscFree(gpiece));
    858:   PetscCall(PetscFree(buffer));
    859:   PetscCall(PetscFPrintf(comm, fp, "\n  </AppendedData>\n"));
    860:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
    861:   PetscCall(PetscFClose(comm, fp));
    862: finalize:
    863:   /* this code sends to rank 0 that writes.
    864:      It may lead to very unbalanced log_view timings
    865:      of the next PETSc function logged.
    866:      Since this call is not performance critical, we
    867:      issue a barrier here to synchronize the processes */
    868:   PetscCall(PetscBarrier((PetscObject)viewer));
    869:   PetscFunctionReturn(PETSC_SUCCESS);
    870: }