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 = PetscDefined(USE_COMPLEX) ? 2 : 1;

    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 (comm == MPI_COMM_NULL) goto finalize;
    195:   PetscCallMPI(MPI_Comm_size(comm, &size));
    196:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

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

    203:   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
    204:   piece.nvertices = 0;
    205:   piece.ncells    = 0;
    206:   piece.nconn     = 0;
    207:   if (!localized) piece.nvertices = vEnd - vStart;
    208:   for (c = cStart; c < cEnd; ++c) {
    209:     PetscInt dof = 0;
    210:     if (hasLabel) {
    211:       PetscInt value;

    213:       PetscCall(DMGetLabelValue(dm, "vtk", c, &value));
    214:       if (value != 1) continue;
    215:     }
    216:     if (localized) PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
    217:     if (!dof) {
    218:       PetscInt *closure = NULL;
    219:       PetscInt  closureSize;

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

    238:   /*
    239:    * Write file header
    240:    */
    241:   if (rank == 0) {
    242:     PetscInt64 boffset = 0;

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

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

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

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

    430:   PetscCall(PetscFPrintf(comm, fp, "  </UnstructuredGrid>\n"));
    431:   PetscCall(PetscFPrintf(comm, fp, "  <AppendedData encoding=\"raw\">\n"));
    432:   PetscCall(PetscFPrintf(comm, fp, "_"));

    434:   if (rank == 0) {
    435:     PetscInt maxsize = 0;
    436:     for (PetscMPIInt r = 0; r < size; r++) {
    437:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nvertices * 3 * sizeof(PetscVTUReal)));
    438:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].ncells * 3 * sizeof(PetscVTUReal)));
    439:       maxsize = PetscMax(maxsize, (PetscInt)(gpiece[r].nconn * sizeof(PetscVTKInt)));
    440:     }
    441:     PetscCall(PetscMalloc(maxsize, &buffer));
    442:   }
    443:   for (PetscMPIInt r = 0; r < size; r++) {
    444:     if (r == rank) {
    445:       PetscInt nsend;
    446:       { /* Position */
    447:         const PetscScalar *x, *cx = NULL;
    448:         PetscVTUReal      *y = NULL;
    449:         Vec                coords, cellCoords;
    450:         const PetscBool    copy = PetscDefined(USE_COMPLEX) ? PETSC_TRUE : (PetscBool)(dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));

    452:         PetscCall(DMGetCoordinatesLocal(dm, &coords));
    453:         PetscCall(VecGetArrayRead(coords, &x));
    454:         PetscCall(DMGetCellCoordinatesLocal(dm, &cellCoords));
    455:         if (cellCoords) PetscCall(VecGetArrayRead(cellCoords, &cx));
    456:         if (copy) {
    457:           PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
    458:           if (localized) {
    459:             PetscInt cnt;
    460:             for (c = cStart, cnt = 0; c < cEnd; c++) {
    461:               PetscInt off, dof;

    463:               PetscCall(PetscSectionGetDof(cellCoordSection, c, &dof));
    464:               if (!dof) {
    465:                 PetscInt *closure = NULL;
    466:                 PetscInt  closureSize;

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

    528:         PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &orank));
    529:         PetscCall(PetscMalloc1(piece.ncells, &owners));
    530:         for (i = 0; i < piece.ncells; i++) owners[i] = orank;
    531:         PetscCall(TransferWrite(comm, viewer, fp, r, 0, owners, buffer, piece.ncells, MPI_INT, tag));
    532:         PetscCall(PetscFree(owners));
    533:       }
    534:       /* Cell data */
    535:       for (link = vtk->link; link; link = link->next) {
    536:         Vec                X = (Vec)link->vec;
    537:         DM                 dmX;
    538:         const PetscScalar *x;
    539:         PetscVTUReal      *y;
    540:         PetscInt           bs      = 1, nfields, field;
    541:         PetscSection       section = NULL;

    543:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
    544:         PetscCall(VecGetDM(X, &dmX));
    545:         if (!dmX) dmX = dm;
    546:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    547:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    548:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
    549:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    550:         field = 0;
    551:         if (link->field >= 0) {
    552:           field   = link->field;
    553:           nfields = field + 1;
    554:         }
    555:         PetscCall(VecGetArrayRead(X, &x));
    556:         PetscCall(PetscMalloc1(piece.ncells * 3, &y));
    557:         for (; field < (nfields ? nfields : 1); field++) {
    558:           PetscInt  fbs, j;
    559:           PetscFV   fv = NULL;
    560:           PetscBool vector;

    562:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
    563:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
    564:           } else fbs = bs; /* Say we have one field with 'bs' components */
    565:           PetscCall(DMGetFieldIfFV_Private(dmX, field, &fv));
    566:           vector = PETSC_FALSE;
    567:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
    568:             vector = PETSC_TRUE;
    569:             for (j = 0; j < fbs; j++) {
    570:               const char *compName = NULL;
    571:               if (fv) {
    572:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
    573:                 if (compName) break;
    574:               }
    575:             }
    576:             if (j < fbs) vector = PETSC_FALSE;
    577:           }
    578:           if (vector) {
    579:             PetscInt cnt, l;
    580:             for (l = 0; l < loops_per_scalar; l++) {
    581:               for (c = cStart, cnt = 0; c < cEnd; c++) {
    582:                 const PetscScalar *xpoint;
    583:                 PetscInt           off, j;

    585:                 if (hasLabel) { /* Ignore some cells */
    586:                   PetscInt value;
    587:                   PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
    588:                   if (value != 1) continue;
    589:                 }
    590:                 if (nfields) PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
    591:                 else PetscCall(PetscSectionGetOffset(section, c, &off));
    592:                 xpoint = &x[off];
    593:                 for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
    594:                 for (; j < 3; j++) y[cnt++] = 0.;
    595:               }
    596:               PetscCheck(cnt == piece.ncells * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    597:               PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells * 3, MPIU_VTUREAL, tag));
    598:             }
    599:           } else {
    600:             for (i = 0; i < fbs; i++) {
    601:               PetscInt cnt, l;
    602:               for (l = 0; l < loops_per_scalar; l++) {
    603:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
    604:                   const PetscScalar *xpoint;
    605:                   PetscInt           off;

    607:                   if (hasLabel) { /* Ignore some cells */
    608:                     PetscInt value;
    609:                     PetscCall(DMGetLabelValue(dmX, "vtk", c, &value));
    610:                     if (value != 1) continue;
    611:                   }
    612:                   if (nfields) PetscCall(PetscSectionGetFieldOffset(section, c, field, &off));
    613:                   else PetscCall(PetscSectionGetOffset(section, c, &off));
    614:                   xpoint   = &x[off];
    615:                   y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
    616:                 }
    617:                 PetscCheck(cnt == piece.ncells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    618:                 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.ncells, MPIU_VTUREAL, tag));
    619:               }
    620:             }
    621:           }
    622:         }
    623:         PetscCall(PetscFree(y));
    624:         PetscCall(VecRestoreArrayRead(X, &x));
    625:       }
    626:       /* point data */
    627:       for (link = vtk->link; link; link = link->next) {
    628:         Vec                X = (Vec)link->vec;
    629:         DM                 dmX;
    630:         const PetscScalar *x;
    631:         PetscVTUReal      *y;
    632:         PetscInt           bs      = 1, nfields, field;
    633:         PetscSection       section = NULL;

    635:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
    636:         PetscCall(VecGetDM(X, &dmX));
    637:         if (!dmX) dmX = dm;
    638:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    639:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    640:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
    641:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    642:         field = 0;
    643:         if (link->field >= 0) {
    644:           field   = link->field;
    645:           nfields = field + 1;
    646:         }
    647:         PetscCall(VecGetArrayRead(X, &x));
    648:         PetscCall(PetscMalloc1(piece.nvertices * 3, &y));
    649:         for (; field < (nfields ? nfields : 1); field++) {
    650:           PetscInt fbs, j;
    651:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
    652:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
    653:           } else fbs = bs; /* Say we have one field with 'bs' components */
    654:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
    655:             PetscInt cnt, l;
    656:             for (l = 0; l < loops_per_scalar; l++) {
    657:               if (!localized) {
    658:                 for (v = vStart, cnt = 0; v < vEnd; v++) {
    659:                   PetscInt           off;
    660:                   const PetscScalar *xpoint;

    662:                   if (nfields) PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
    663:                   else PetscCall(PetscSectionGetOffset(section, v, &off));
    664:                   xpoint = &x[off];
    665:                   for (j = 0; j < fbs; j++) y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
    666:                   for (; j < 3; j++) y[cnt++] = 0.;
    667:                 }
    668:               } else {
    669:                 for (c = cStart, cnt = 0; c < cEnd; c++) {
    670:                   PetscInt *closure = NULL;
    671:                   PetscInt  closureSize, off;

    673:                   PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    674:                   for (v = 0, off = 0; v < closureSize * 2; v += 2) {
    675:                     if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
    676:                       PetscInt           voff;
    677:                       const PetscScalar *xpoint;

    679:                       if (nfields) {
    680:                         PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
    681:                       } else {
    682:                         PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
    683:                       }
    684:                       xpoint = &x[voff];
    685:                       for (j = 0; j < fbs; j++) y[cnt + off++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
    686:                       for (; j < 3; j++) y[cnt + off++] = 0.;
    687:                     }
    688:                   }
    689:                   cnt += off;
    690:                   PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    691:                 }
    692:               }
    693:               PetscCheck(cnt == piece.nvertices * 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    694:               PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices * 3, MPIU_VTUREAL, tag));
    695:             }
    696:           } else {
    697:             for (i = 0; i < fbs; i++) {
    698:               PetscInt cnt, l;
    699:               for (l = 0; l < loops_per_scalar; l++) {
    700:                 if (!localized) {
    701:                   for (v = vStart, cnt = 0; v < vEnd; v++) {
    702:                     PetscInt           off;
    703:                     const PetscScalar *xpoint;

    705:                     if (nfields) PetscCall(PetscSectionGetFieldOffset(section, v, field, &off));
    706:                     else PetscCall(PetscSectionGetOffset(section, v, &off));
    707:                     xpoint   = &x[off];
    708:                     y[cnt++] = (PetscVTUReal)(l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
    709:                   }
    710:                 } else {
    711:                   for (c = cStart, cnt = 0; c < cEnd; c++) {
    712:                     PetscInt *closure = NULL;
    713:                     PetscInt  closureSize, off;

    715:                     PetscCall(DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    716:                     for (v = 0, off = 0; v < closureSize * 2; v += 2) {
    717:                       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
    718:                         PetscInt           voff;
    719:                         const PetscScalar *xpoint;

    721:                         if (nfields) PetscCall(PetscSectionGetFieldOffset(section, closure[v], field, &voff));
    722:                         else PetscCall(PetscSectionGetOffset(section, closure[v], &voff));
    723:                         xpoint         = &x[voff];
    724:                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
    725:                       }
    726:                     }
    727:                     cnt += off;
    728:                     PetscCall(DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure));
    729:                   }
    730:                 }
    731:                 PetscCheck(cnt == piece.nvertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count does not match");
    732:                 PetscCall(TransferWrite(comm, viewer, fp, r, 0, y, buffer, piece.nvertices, MPIU_VTUREAL, tag));
    733:               }
    734:             }
    735:           }
    736:         }
    737:         PetscCall(PetscFree(y));
    738:         PetscCall(VecRestoreArrayRead(X, &x));
    739:       }
    740:     } else if (rank == 0) {
    741:       PetscInt l;

    743:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag)); /* positions */
    744:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nconn, MPI_INT, tag));              /* connectivity */
    745:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag));             /* offsets */
    746:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_CHAR, tag));            /* types */
    747:       PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPI_INT, tag));             /* owner rank (cells) */
    748:       /* all cell data */
    749:       for (link = vtk->link; link; link = link->next) {
    750:         Vec          X  = (Vec)link->vec;
    751:         PetscInt     bs = 1, nfields, field;
    752:         DM           dmX;
    753:         PetscSection section = NULL;

    755:         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
    756:         PetscCall(VecGetDM(X, &dmX));
    757:         if (!dmX) dmX = dm;
    758:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    759:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    760:         if (cEnd > cStart) PetscCall(PetscSectionGetDof(section, cStart, &bs));
    761:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    762:         field = 0;
    763:         if (link->field >= 0) {
    764:           field   = link->field;
    765:           nfields = field + 1;
    766:         }
    767:         for (; field < (nfields ? nfields : 1); field++) {
    768:           PetscInt  fbs, j;
    769:           PetscFV   fv = NULL;
    770:           PetscBool vector;

    772:           if (nfields && cEnd > cStart) { /* We have user-defined fields/components */
    773:             PetscCall(PetscSectionGetFieldDof(section, cStart, field, &fbs));
    774:           } else fbs = bs; /* Say we have one field with 'bs' components */
    775:           PetscCall(DMGetFieldIfFV_Private(dmX, field, &fv));
    776:           vector = PETSC_FALSE;
    777:           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
    778:             vector = PETSC_TRUE;
    779:             for (j = 0; j < fbs; j++) {
    780:               const char *compName = NULL;
    781:               if (fv) {
    782:                 PetscCall(PetscFVGetComponentName(fv, j, &compName));
    783:                 if (compName) break;
    784:               }
    785:             }
    786:             if (j < fbs) vector = PETSC_FALSE;
    787:           }
    788:           if (vector) {
    789:             for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells * 3, MPIU_VTUREAL, tag));
    790:           } else {
    791:             for (i = 0; i < fbs; i++) {
    792:               for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].ncells, MPIU_VTUREAL, tag));
    793:             }
    794:           }
    795:         }
    796:       }
    797:       /* all point data */
    798:       for (link = vtk->link; link; link = link->next) {
    799:         Vec          X = (Vec)link->vec;
    800:         DM           dmX;
    801:         PetscInt     bs      = 1, nfields, field;
    802:         PetscSection section = NULL;

    804:         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
    805:         PetscCall(VecGetDM(X, &dmX));
    806:         if (!dmX) dmX = dm;
    807:         PetscCall(PetscObjectQuery(link->vec, "section", (PetscObject *)&section));
    808:         if (!section) PetscCall(DMGetLocalSection(dmX, &section));
    809:         if (vEnd > vStart) PetscCall(PetscSectionGetDof(section, vStart, &bs));
    810:         PetscCall(PetscSectionGetNumFields(section, &nfields));
    811:         field = 0;
    812:         if (link->field >= 0) {
    813:           field   = link->field;
    814:           nfields = field + 1;
    815:         }
    816:         for (; field < (nfields ? nfields : 1); field++) {
    817:           PetscInt fbs;
    818:           if (nfields && vEnd > vStart) { /* We have user-defined fields/components */
    819:             PetscCall(PetscSectionGetFieldDof(section, vStart, field, &fbs));
    820:           } else fbs = bs; /* Say we have one field with 'bs' components */
    821:           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
    822:             for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices * 3, MPIU_VTUREAL, tag));
    823:           } else {
    824:             for (i = 0; i < fbs; i++) {
    825:               for (l = 0; l < loops_per_scalar; l++) PetscCall(TransferWrite(comm, viewer, fp, r, 0, NULL, buffer, gpiece[r].nvertices, MPIU_VTUREAL, tag));
    826:             }
    827:           }
    828:         }
    829:       }
    830:     }
    831:   }
    832:   PetscCall(PetscFree(gpiece));
    833:   PetscCall(PetscFree(buffer));
    834:   PetscCall(PetscFPrintf(comm, fp, "\n  </AppendedData>\n"));
    835:   PetscCall(PetscFPrintf(comm, fp, "</VTKFile>\n"));
    836:   PetscCall(PetscFClose(comm, fp));
    837: finalize:
    838:   /* this code sends to rank 0 that writes.
    839:      It may lead to very unbalanced log_view timings
    840:      of the next PETSc function logged.
    841:      Since this call is not performance critical, we
    842:      issue a barrier here to synchronize the processes */
    843:   PetscCall(PetscBarrier((PetscObject)viewer));
    844:   PetscFunctionReturn(PETSC_SUCCESS);
    845: }