Actual source code: plexhdf5.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/isimpl.h>
      3: #include <petsc/private/vecimpl.h>
      4: #include <petsc/private/viewerhdf5impl.h>
      5: #include <petsclayouthdf5.h>

      7: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

      9: #if defined(PETSC_HAVE_HDF5)
     10: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
     11: {
     12:   Vec            stamp;
     13:   PetscMPIInt    rank;

     17:   if (seqnum < 0) return(0);
     18:   MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
     19:   VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
     20:   VecSetBlockSize(stamp, 1);
     21:   PetscObjectSetName((PetscObject) stamp, seqname);
     22:   if (!rank) {
     23:     PetscReal timeScale;
     24:     PetscBool istime;

     26:     PetscStrncmp(seqname, "time", 5, &istime);
     27:     if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); value *= timeScale;}
     28:     VecSetValue(stamp, 0, value, INSERT_VALUES);
     29:   }
     30:   VecAssemblyBegin(stamp);
     31:   VecAssemblyEnd(stamp);
     32:   PetscViewerHDF5PushGroup(viewer, "/");
     33:   PetscViewerHDF5SetTimestep(viewer, seqnum);
     34:   VecView(stamp, viewer);
     35:   PetscViewerHDF5PopGroup(viewer);
     36:   VecDestroy(&stamp);
     37:   return(0);
     38: }

     40: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
     41: {
     42:   Vec            stamp;
     43:   PetscMPIInt    rank;

     47:   if (seqnum < 0) return(0);
     48:   MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
     49:   VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
     50:   VecSetBlockSize(stamp, 1);
     51:   PetscObjectSetName((PetscObject) stamp, seqname);
     52:   PetscViewerHDF5PushGroup(viewer, "/");
     53:   PetscViewerHDF5SetTimestep(viewer, seqnum);
     54:   VecLoad(stamp, viewer);
     55:   PetscViewerHDF5PopGroup(viewer);
     56:   if (!rank) {
     57:     const PetscScalar *a;
     58:     PetscReal timeScale;
     59:     PetscBool istime;

     61:     VecGetArrayRead(stamp, &a);
     62:     *value = a[0];
     63:     VecRestoreArrayRead(stamp, &a);
     64:     PetscStrncmp(seqname, "time", 5, &istime);
     65:     if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
     66:   }
     67:   VecDestroy(&stamp);
     68:   return(0);
     69: }

     71: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
     72: {
     73:   IS              cutcells = NULL;
     74:   const PetscInt *cutc;
     75:   PetscInt        cellHeight, vStart, vEnd, cStart, cEnd, c;
     76:   PetscErrorCode  ierr;

     79:   if (!cutLabel) return(0);
     80:   DMPlexGetVTKCellHeight(dm, &cellHeight);
     81:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
     82:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
     83:   /* Label vertices that should be duplicated */
     84:   DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
     85:   DMLabelGetStratumIS(cutLabel, 2, &cutcells);
     86:   if (cutcells) {
     87:     PetscInt n;

     89:     ISGetIndices(cutcells, &cutc);
     90:     ISGetLocalSize(cutcells, &n);
     91:     for (c = 0; c < n; ++c) {
     92:       if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
     93:         PetscInt *closure = NULL;
     94:         PetscInt  closureSize, cl, value;

     96:         DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
     97:         for (cl = 0; cl < closureSize*2; cl += 2) {
     98:           if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
     99:             DMLabelGetValue(cutLabel, closure[cl], &value);
    100:             if (value == 1) {
    101:               DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
    102:             }
    103:           }
    104:         }
    105:         DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
    106:       }
    107:     }
    108:     ISRestoreIndices(cutcells, &cutc);
    109:   }
    110:   ISDestroy(&cutcells);
    111:   return(0);
    112: }

    114: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
    115: {
    116:   DM                      dm;
    117:   DM                      dmBC;
    118:   PetscSection            section, sectionGlobal;
    119:   Vec                     gv;
    120:   const char             *name;
    121:   PetscViewerVTKFieldType ft;
    122:   PetscViewerFormat       format;
    123:   PetscInt                seqnum;
    124:   PetscReal               seqval;
    125:   PetscBool               isseq;
    126:   PetscErrorCode          ierr;

    129:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
    130:   VecGetDM(v, &dm);
    131:   DMGetLocalSection(dm, &section);
    132:   DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
    133:   PetscViewerHDF5SetTimestep(viewer, seqnum);
    134:   DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
    135:   PetscViewerGetFormat(viewer, &format);
    136:   DMGetOutputDM(dm, &dmBC);
    137:   DMGetGlobalSection(dmBC, &sectionGlobal);
    138:   DMGetGlobalVector(dmBC, &gv);
    139:   PetscObjectGetName((PetscObject) v, &name);
    140:   PetscObjectSetName((PetscObject) gv, name);
    141:   DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
    142:   DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
    143:   PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
    144:   if (isseq) {VecView_Seq(gv, viewer);}
    145:   else       {VecView_MPI(gv, viewer);}
    146:   if (format == PETSC_VIEWER_HDF5_VIZ) {
    147:     /* Output visualization representation */
    148:     PetscInt numFields, f;
    149:     DMLabel  cutLabel, cutVertexLabel = NULL;

    151:     PetscSectionGetNumFields(section, &numFields);
    152:     DMGetLabel(dm, "periodic_cut", &cutLabel);
    153:     for (f = 0; f < numFields; ++f) {
    154:       Vec         subv;
    155:       IS          is;
    156:       const char *fname, *fgroup;
    157:       char        subname[PETSC_MAX_PATH_LEN];
    158:       PetscInt    pStart, pEnd;

    160:       DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
    161:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
    162:       PetscSectionGetFieldName(section, f, &fname);
    163:       if (!fname) continue;
    164:       PetscViewerHDF5PushGroup(viewer, fgroup);
    165:       if (cutLabel) {
    166:         const PetscScalar *ga;
    167:         PetscScalar       *suba;
    168:         PetscInt           Nc, gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;

    170:         DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
    171:         PetscSectionGetFieldComponents(section, f, &Nc);
    172:         for (p = pStart; p < pEnd; ++p) {
    173:           PetscInt gdof, fdof = 0, val;

    175:           PetscSectionGetDof(sectionGlobal, p, &gdof);
    176:           if (gdof > 0) {PetscSectionGetFieldDof(section, p, f, &fdof);}
    177:           subSize += fdof;
    178:           DMLabelGetValue(cutVertexLabel, p, &val);
    179:           if (val == 1) extSize += fdof;
    180:         }
    181:         VecCreate(PetscObjectComm((PetscObject) gv), &subv);
    182:         VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
    183:         VecSetBlockSize(subv, Nc);
    184:         VecSetType(subv, VECSTANDARD);
    185:         VecGetOwnershipRange(gv, &gstart, NULL);
    186:         VecGetArrayRead(gv, &ga);
    187:         VecGetArray(subv, &suba);
    188:         for (p = pStart; p < pEnd; ++p) {
    189:           PetscInt gdof, goff, val;

    191:           PetscSectionGetDof(sectionGlobal, p, &gdof);
    192:           if (gdof > 0) {
    193:             PetscInt fdof, fc, f2, poff = 0;

    195:             PetscSectionGetOffset(sectionGlobal, p, &goff);
    196:             /* Can get rid of this loop by storing field information in the global section */
    197:             for (f2 = 0; f2 < f; ++f2) {
    198:               PetscSectionGetFieldDof(section, p, f2, &fdof);
    199:               poff += fdof;
    200:             }
    201:             PetscSectionGetFieldDof(section, p, f, &fdof);
    202:             for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
    203:             DMLabelGetValue(cutVertexLabel, p, &val);
    204:             if (val == 1) {
    205:               for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
    206:             }
    207:           }
    208:         }
    209:         VecRestoreArrayRead(gv, &ga);
    210:         VecRestoreArray(subv, &suba);
    211:         DMLabelDestroy(&cutVertexLabel);
    212:       } else {
    213:         PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
    214:       }
    215:       PetscStrncpy(subname, name,sizeof(subname));
    216:       PetscStrlcat(subname, "_",sizeof(subname));
    217:       PetscStrlcat(subname, fname,sizeof(subname));
    218:       PetscObjectSetName((PetscObject) subv, subname);
    219:       if (isseq) {VecView_Seq(subv, viewer);}
    220:       else       {VecView_MPI(subv, viewer);}
    221:       if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
    222:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "vector");
    223:       } else {
    224:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "scalar");
    225:       }
    226:       if (cutLabel) {VecDestroy(&subv);}
    227:       else          {PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);}
    228:       PetscViewerHDF5PopGroup(viewer);
    229:     }
    230:   }
    231:   DMRestoreGlobalVector(dmBC, &gv);
    232:   return(0);
    233: }

    235: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
    236: {
    237:   DM             dm;
    238:   Vec            locv;
    239:   PetscObject    isZero;
    240:   const char    *name;
    241:   PetscReal      time;

    245:   VecGetDM(v, &dm);
    246:   DMGetLocalVector(dm, &locv);
    247:   PetscObjectGetName((PetscObject) v, &name);
    248:   PetscObjectSetName((PetscObject) locv, name);
    249:   PetscObjectQuery((PetscObject) v, "__Vec_bc_zero__", &isZero);
    250:   PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", isZero);
    251:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
    252:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
    253:   DMGetOutputSequenceNumber(dm, NULL, &time);
    254:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
    255:   PetscViewerHDF5PushGroup(viewer, "/fields");
    256:   PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
    257:   VecView_Plex_Local_HDF5_Internal(locv, viewer);
    258:   PetscViewerPopFormat(viewer);
    259:   PetscViewerHDF5PopGroup(viewer);
    260:   PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", NULL);
    261:   DMRestoreLocalVector(dm, &locv);
    262:   return(0);
    263: }

    265: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
    266: {
    267:   PetscBool      isseq;

    271:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
    272:   PetscViewerHDF5PushGroup(viewer, "/fields");
    273:   if (isseq) {VecView_Seq(v, viewer);}
    274:   else       {VecView_MPI(v, viewer);}
    275:   PetscViewerHDF5PopGroup(viewer);
    276:   return(0);
    277: }

    279: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
    280: {
    281:   DM             dm;
    282:   Vec            locv;
    283:   const char    *name;
    284:   PetscInt       seqnum;

    288:   VecGetDM(v, &dm);
    289:   DMGetLocalVector(dm, &locv);
    290:   PetscObjectGetName((PetscObject) v, &name);
    291:   PetscObjectSetName((PetscObject) locv, name);
    292:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
    293:   PetscViewerHDF5SetTimestep(viewer, seqnum);
    294:   PetscViewerHDF5PushGroup(viewer, "/fields");
    295:   VecLoad_Plex_Local(locv, viewer);
    296:   PetscViewerHDF5PopGroup(viewer);
    297:   DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
    298:   DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
    299:   DMRestoreLocalVector(dm, &locv);
    300:   return(0);
    301: }

    303: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
    304: {
    305:   DM             dm;
    306:   PetscInt       seqnum;

    310:   VecGetDM(v, &dm);
    311:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
    312:   PetscViewerHDF5SetTimestep(viewer, seqnum);
    313:   PetscViewerHDF5PushGroup(viewer, "/fields");
    314:   VecLoad_Default(v, viewer);
    315:   PetscViewerHDF5PopGroup(viewer);
    316:   return(0);
    317: }

    319: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
    320: {
    321:   IS              orderIS, conesIS, cellsIS, orntsIS;
    322:   const PetscInt *gpoint;
    323:   PetscInt       *order, *sizes, *cones, *ornts;
    324:   PetscInt        dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
    325:   PetscErrorCode  ierr;

    328:   ISGetIndices(globalPointNumbers, &gpoint);
    329:   DMGetDimension(dm, &dim);
    330:   DMPlexGetChart(dm, &pStart, &pEnd);
    331:   for (p = pStart; p < pEnd; ++p) {
    332:     if (gpoint[p] >= 0) {
    333:       PetscInt coneSize;

    335:       DMPlexGetConeSize(dm, p, &coneSize);
    336:       conesSize += 1;
    337:       cellsSize += coneSize;
    338:     }
    339:   }
    340:   PetscMalloc1(conesSize, &order);
    341:   PetscMalloc1(conesSize, &sizes);
    342:   PetscMalloc1(cellsSize, &cones);
    343:   PetscMalloc1(cellsSize, &ornts);
    344:   for (p = pStart; p < pEnd; ++p) {
    345:     if (gpoint[p] >= 0) {
    346:       const PetscInt *cone, *ornt;
    347:       PetscInt        coneSize, cp;

    349:       DMPlexGetConeSize(dm, p, &coneSize);
    350:       DMPlexGetCone(dm, p, &cone);
    351:       DMPlexGetConeOrientation(dm, p, &ornt);
    352:       order[s]   = gpoint[p];
    353:       sizes[s++] = coneSize;
    354:       for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
    355:     }
    356:   }
    357:   if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %D != %D", s, conesSize);
    358:   if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %D != %D", c, cellsSize);
    359:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
    360:   PetscObjectSetName((PetscObject) orderIS, "order");
    361:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
    362:   PetscObjectSetName((PetscObject) conesIS, "cones");
    363:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
    364:   PetscObjectSetName((PetscObject) cellsIS, "cells");
    365:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
    366:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
    367:   PetscViewerHDF5PushGroup(viewer, "/topology");
    368:   ISView(orderIS, viewer);
    369:   ISView(conesIS, viewer);
    370:   ISView(cellsIS, viewer);
    371:   PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
    372:   ISView(orntsIS, viewer);
    373:   PetscViewerHDF5PopGroup(viewer);
    374:   ISDestroy(&orderIS);
    375:   ISDestroy(&conesIS);
    376:   ISDestroy(&cellsIS);
    377:   ISDestroy(&orntsIS);
    378:   ISRestoreIndices(globalPointNumbers, &gpoint);
    379:   return(0);
    380: }

    382: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
    383: {
    384:   PetscSF         sfPoint;
    385:   DMLabel         cutLabel, cutVertexLabel = NULL;
    386:   IS              globalVertexNumbers, cutvertices = NULL;
    387:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
    388:   PetscInt       *vertices;
    389:   PetscInt        conesSize = 0;
    390:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
    391:   PetscErrorCode  ierr;

    394:   *numCorners = 0;
    395:   DMGetDimension(dm, &dim);
    396:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    397:   ISGetIndices(globalCellNumbers, &gcell);

    399:   for (cell = cStart; cell < cEnd; ++cell) {
    400:     PetscInt *closure = NULL;
    401:     PetscInt  closureSize, v, Nc = 0;

    403:     if (gcell[cell] < 0) continue;
    404:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    405:     for (v = 0; v < closureSize*2; v += 2) {
    406:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
    407:     }
    408:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    409:     conesSize += Nc;
    410:     if (!numCornersLocal)           numCornersLocal = Nc;
    411:     else if (numCornersLocal != Nc) numCornersLocal = 1;
    412:   }
    413:   MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
    414:   if (numCornersLocal && (numCornersLocal != *numCorners || *numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
    415:   /* Handle periodic cuts by identifying vertices which should be duplicated */
    416:   DMGetLabel(dm, "periodic_cut", &cutLabel);
    417:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
    418:   if (cutVertexLabel) {DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);}
    419:   if (cutvertices) {
    420:     ISGetIndices(cutvertices, &cutverts);
    421:     ISGetLocalSize(cutvertices, &vExtra);
    422:   }
    423:   DMGetPointSF(dm, &sfPoint);
    424:   if (cutLabel) {
    425:     const PetscInt    *ilocal;
    426:     const PetscSFNode *iremote;
    427:     PetscInt           nroots, nleaves;

    429:     PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
    430:     if (nleaves < 0) {
    431:       PetscObjectReference((PetscObject) sfPoint);
    432:     } else {
    433:       PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
    434:       PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, ilocal, PETSC_USE_POINTER, iremote, PETSC_USE_POINTER);
    435:     }
    436:   } else {
    437:     PetscObjectReference((PetscObject) sfPoint);
    438:   }
    439:   /* Number all vertices */
    440:   DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
    441:   PetscSFDestroy(&sfPoint);
    442:   /* Create cones */
    443:   ISGetIndices(globalVertexNumbers, &gvertex);
    444:   PetscMalloc1(conesSize, &vertices);
    445:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
    446:     PetscInt *closure = NULL;
    447:     PetscInt  closureSize, Nc = 0, p, value = -1;
    448:     PetscBool replace;

    450:     if (gcell[cell] < 0) continue;
    451:     if (cutLabel) {DMLabelGetValue(cutLabel, cell, &value);}
    452:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
    453:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    454:     for (p = 0; p < closureSize*2; p += 2) {
    455:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
    456:         closure[Nc++] = closure[p];
    457:       }
    458:     }
    459:     DMPlexReorderCell(dm, cell, closure);
    460:     for (p = 0; p < Nc; ++p) {
    461:       PetscInt nv, gv = gvertex[closure[p] - vStart];

    463:       if (replace) {
    464:         PetscFindInt(closure[p], vExtra, cutverts, &nv);
    465:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
    466:       }
    467:       vertices[v++] = gv < 0 ? -(gv+1) : gv;
    468:     }
    469:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
    470:   }
    471:   ISRestoreIndices(globalVertexNumbers, &gvertex);
    472:   ISDestroy(&globalVertexNumbers);
    473:   ISRestoreIndices(globalCellNumbers, &gcell);
    474:   if (cutvertices) {ISRestoreIndices(cutvertices, &cutverts);}
    475:   ISDestroy(&cutvertices);
    476:   DMLabelDestroy(&cutVertexLabel);
    477:   if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %D != %D", v, conesSize);
    478:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
    479:   PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
    480:   PetscObjectSetName((PetscObject) *cellIS, "cells");
    481:   return(0);
    482: }

    484: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
    485: {
    486:   DM              cdm;
    487:   DMLabel         depthLabel, ctLabel;
    488:   IS              cellIS;
    489:   PetscInt        dim, depth, cellHeight, c;
    490:   hid_t           fileId, groupId;
    491:   PetscErrorCode  ierr;

    494:   PetscViewerHDF5PushGroup(viewer, "/viz");
    495:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
    496:   PetscStackCallHDF5(H5Gclose,(groupId));

    498:   PetscViewerHDF5PopGroup(viewer);
    499:   DMGetDimension(dm, &dim);
    500:   DMPlexGetDepth(dm, &depth);
    501:   DMGetCoordinateDM(dm, &cdm);
    502:   DMPlexGetVTKCellHeight(dm, &cellHeight);
    503:   DMPlexGetDepthLabel(dm, &depthLabel);
    504:   DMPlexGetCellTypeLabel(dm, &ctLabel);
    505:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    506:     const DMPolytopeType ict = (DMPolytopeType) c;
    507:     PetscInt             pStart, pEnd, dep, numCorners, n = 0;
    508:     PetscBool            output = PETSC_FALSE, doOutput;

    510:     if (ict == DM_POLYTOPE_FV_GHOST) continue;
    511:     DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
    512:     if (pStart >= 0) {
    513:       DMLabelGetValue(depthLabel, pStart, &dep);
    514:       if (dep == depth - cellHeight) output = PETSC_TRUE;
    515:     }
    516:     MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
    517:     if (!doOutput) continue;
    518:     CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners,  &cellIS);
    519:     if (!n) {
    520:       PetscViewerHDF5PushGroup(viewer, "/viz/topology");
    521:     } else {
    522:       char group[PETSC_MAX_PATH_LEN];

    524:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
    525:       PetscViewerHDF5PushGroup(viewer, group);
    526:     }
    527:     ISView(cellIS, viewer);
    528:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
    529:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim",     PETSC_INT, (void *) &dim);
    530:     ISDestroy(&cellIS);
    531:     PetscViewerHDF5PopGroup(viewer);
    532:     ++n;
    533:   }
    534:   return(0);
    535: }

    537: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
    538: {
    539:   DM             cdm;
    540:   Vec            coordinates, newcoords;
    541:   PetscReal      lengthScale;
    542:   PetscInt       m, M, bs;

    546:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
    547:   DMGetCoordinateDM(dm, &cdm);
    548:   DMGetCoordinates(dm, &coordinates);
    549:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
    550:   PetscObjectSetName((PetscObject) newcoords, "vertices");
    551:   VecGetSize(coordinates, &M);
    552:   VecGetLocalSize(coordinates, &m);
    553:   VecSetSizes(newcoords, m, M);
    554:   VecGetBlockSize(coordinates, &bs);
    555:   VecSetBlockSize(newcoords, bs);
    556:   VecSetType(newcoords,VECSTANDARD);
    557:   VecCopy(coordinates, newcoords);
    558:   VecScale(newcoords, lengthScale);
    559:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
    560:   PetscViewerHDF5PushGroup(viewer, "/geometry");
    561:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
    562:   VecView(newcoords, viewer);
    563:   PetscViewerPopFormat(viewer);
    564:   PetscViewerHDF5PopGroup(viewer);
    565:   VecDestroy(&newcoords);
    566:   return(0);
    567: }

    569: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
    570: {
    571:   DM               cdm;
    572:   Vec              coordinatesLocal, newcoords;
    573:   PetscSection     cSection, cGlobalSection;
    574:   PetscScalar     *coords, *ncoords;
    575:   DMLabel          cutLabel, cutVertexLabel = NULL;
    576:   const PetscReal *L;
    577:   const DMBoundaryType *bd;
    578:   PetscReal        lengthScale;
    579:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
    580:   PetscBool        localized, embedded;
    581:   hid_t            fileId, groupId;
    582:   PetscErrorCode   ierr;

    585:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    586:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
    587:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
    588:   VecGetBlockSize(coordinatesLocal, &bs);
    589:   DMGetCoordinatesLocalized(dm, &localized);
    590:   if (localized == PETSC_FALSE) return(0);
    591:   DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
    592:   DMGetCoordinateDM(dm, &cdm);
    593:   DMGetLocalSection(cdm, &cSection);
    594:   DMGetGlobalSection(cdm, &cGlobalSection);
    595:   DMGetLabel(dm, "periodic_cut", &cutLabel);
    596:   N    = 0;

    598:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
    599:   VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
    600:   PetscSectionGetDof(cSection, vStart, &dof);
    601:   PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
    602:   embedded  = (PetscBool) (L && dof == 2 && !cutLabel);
    603:   if (cutVertexLabel) {
    604:     DMLabelGetStratumSize(cutVertexLabel, 1, &v);
    605:     N   += dof*v;
    606:   }
    607:   for (v = vStart; v < vEnd; ++v) {
    608:     PetscSectionGetDof(cGlobalSection, v, &dof);
    609:     if (dof < 0) continue;
    610:     if (embedded) N += dof+1;
    611:     else          N += dof;
    612:   }
    613:   if (embedded) {VecSetBlockSize(newcoords, bs+1);}
    614:   else          {VecSetBlockSize(newcoords, bs);}
    615:   VecSetSizes(newcoords, N, PETSC_DETERMINE);
    616:   VecSetType(newcoords, VECSTANDARD);
    617:   VecGetArray(coordinatesLocal, &coords);
    618:   VecGetArray(newcoords,        &ncoords);
    619:   coordSize = 0;
    620:   for (v = vStart; v < vEnd; ++v) {
    621:     PetscSectionGetDof(cGlobalSection, v, &dof);
    622:     PetscSectionGetOffset(cSection, v, &off);
    623:     if (dof < 0) continue;
    624:     if (embedded) {
    625:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
    626:         PetscReal theta, phi, r, R;
    627:         /* XY-periodic */
    628:         /* Suppose its an y-z circle, then
    629:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
    630:            and the circle in that plane is
    631:              \hat r cos(phi) + \hat x sin(phi) */
    632:         theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
    633:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
    634:         r     = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
    635:         R     = L[1]/(2.0*PETSC_PI);
    636:         ncoords[coordSize++] =  PetscSinReal(phi) * r;
    637:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
    638:         ncoords[coordSize++] =  PetscSinReal(theta) * (R + r * PetscCosReal(phi));
    639:       } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
    640:         /* X-periodic */
    641:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
    642:         ncoords[coordSize++] = coords[off+1];
    643:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
    644:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
    645:         /* Y-periodic */
    646:         ncoords[coordSize++] = coords[off+0];
    647:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
    648:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
    649:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
    650:         PetscReal phi, r, R;
    651:         /* Mobius strip */
    652:         /* Suppose its an x-z circle, then
    653:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
    654:            and in that plane we rotate by pi as we go around the circle
    655:              \hat r cos(phi/2) + \hat y sin(phi/2) */
    656:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
    657:         R     = L[0];
    658:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
    659:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
    660:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
    661:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
    662:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
    663:     } else {
    664:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
    665:     }
    666:   }
    667:   if (cutVertexLabel) {
    668:     IS              vertices;
    669:     const PetscInt *verts;
    670:     PetscInt        n;

    672:     DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
    673:     if (vertices) {
    674:       ISGetIndices(vertices, &verts);
    675:       ISGetLocalSize(vertices, &n);
    676:       for (v = 0; v < n; ++v) {
    677:         PetscSectionGetDof(cSection, verts[v], &dof);
    678:         PetscSectionGetOffset(cSection, verts[v], &off);
    679:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
    680:       }
    681:       ISRestoreIndices(vertices, &verts);
    682:       ISDestroy(&vertices);
    683:     }
    684:   }
    685:   if (coordSize != N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %D != %D", coordSize, N);
    686:   DMLabelDestroy(&cutVertexLabel);
    687:   VecRestoreArray(coordinatesLocal, &coords);
    688:   VecRestoreArray(newcoords,        &ncoords);
    689:   PetscObjectSetName((PetscObject) newcoords, "vertices");
    690:   VecScale(newcoords, lengthScale);
    691:   PetscViewerHDF5PushGroup(viewer, "/viz");
    692:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
    693:   PetscStackCallHDF5(H5Gclose,(groupId));
    694:   PetscViewerHDF5PopGroup(viewer);
    695:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
    696:   VecView(newcoords, viewer);
    697:   PetscViewerHDF5PopGroup(viewer);
    698:   VecDestroy(&newcoords);
    699:   return(0);
    700: }

    703: static PetscErrorCode DMPlexWriteLabels_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
    704: {
    705:   const PetscInt   *gpoint;
    706:   PetscInt          numLabels, l;
    707:   hid_t             fileId, groupId;
    708:   PetscErrorCode    ierr;

    711:   ISGetIndices(globalPointNumbers, &gpoint);
    712:   PetscViewerHDF5PushGroup(viewer, "/labels");
    713:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
    714:   if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
    715:   PetscViewerHDF5PopGroup(viewer);
    716:   DMGetNumLabels(dm, &numLabels);
    717:   for (l = 0; l < numLabels; ++l) {
    718:     DMLabel         label;
    719:     const char     *name;
    720:     IS              valueIS, pvalueIS, globalValueIS;
    721:     const PetscInt *values;
    722:     PetscInt        numValues, v;
    723:     PetscBool       isDepth, output;
    724:     char            group[PETSC_MAX_PATH_LEN];

    726:     DMGetLabelName(dm, l, &name);
    727:     DMGetLabelOutput(dm, name, &output);
    728:     PetscStrncmp(name, "depth", 10, &isDepth);
    729:     if (isDepth || !output) continue;
    730:     DMGetLabel(dm, name, &label);
    731:     PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
    732:     PetscViewerHDF5PushGroup(viewer, group);
    733:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
    734:     if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
    735:     PetscViewerHDF5PopGroup(viewer);
    736:     DMLabelGetValueIS(label, &valueIS);
    737:     /* Must copy to a new IS on the global comm */
    738:     ISGetLocalSize(valueIS, &numValues);
    739:     ISGetIndices(valueIS, &values);
    740:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
    741:     ISRestoreIndices(valueIS, &values);
    742:     ISAllGather(pvalueIS, &globalValueIS);
    743:     ISDestroy(&pvalueIS);
    744:     ISSortRemoveDups(globalValueIS);
    745:     ISGetLocalSize(globalValueIS, &numValues);
    746:     ISGetIndices(globalValueIS, &values);
    747:     for (v = 0; v < numValues; ++v) {
    748:       IS              stratumIS, globalStratumIS;
    749:       const PetscInt *spoints = NULL;
    750:       PetscInt       *gspoints, n = 0, gn, p;
    751:       const char     *iname = "indices";

    753:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%D", name, values[v]);
    754:       DMLabelGetStratumIS(label, values[v], &stratumIS);

    756:       if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
    757:       if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
    758:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
    759:       PetscMalloc1(gn,&gspoints);
    760:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
    761:       if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
    762:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
    763:       if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
    764:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

    766:       PetscViewerHDF5PushGroup(viewer, group);
    767:       ISView(globalStratumIS, viewer);
    768:       PetscViewerHDF5PopGroup(viewer);
    769:       ISDestroy(&globalStratumIS);
    770:       ISDestroy(&stratumIS);
    771:     }
    772:     ISRestoreIndices(globalValueIS, &values);
    773:     ISDestroy(&globalValueIS);
    774:     ISDestroy(&valueIS);
    775:   }
    776:   ISRestoreIndices(globalPointNumbers, &gpoint);
    777:   return(0);
    778: }

    780: /* We only write cells and vertices. Does this screw up parallel reading? */
    781: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
    782: {
    783:   IS                globalPointNumbers;
    784:   PetscViewerFormat format;
    785:   PetscBool         viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
    786:   PetscErrorCode    ierr;

    789:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
    790:   DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
    791:   DMPlexWriteLabels_HDF5_Static(dm, globalPointNumbers, viewer);

    793:   PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);

    795:   PetscViewerGetFormat(viewer, &format);
    796:   switch (format) {
    797:     case PETSC_VIEWER_HDF5_VIZ:
    798:       viz_geom    = PETSC_TRUE;
    799:       xdmf_topo   = PETSC_TRUE;
    800:       break;
    801:     case PETSC_VIEWER_HDF5_XDMF:
    802:       xdmf_topo   = PETSC_TRUE;
    803:       break;
    804:     case PETSC_VIEWER_HDF5_PETSC:
    805:       petsc_topo  = PETSC_TRUE;
    806:       break;
    807:     case PETSC_VIEWER_DEFAULT:
    808:     case PETSC_VIEWER_NATIVE:
    809:       viz_geom    = PETSC_TRUE;
    810:       xdmf_topo   = PETSC_TRUE;
    811:       petsc_topo  = PETSC_TRUE;
    812:       break;
    813:     default:
    814:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
    815:   }

    817:   if (viz_geom)   {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
    818:   if (xdmf_topo)  {DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);}
    819:   if (petsc_topo) {DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);}

    821:   ISDestroy(&globalPointNumbers);
    822:   return(0);
    823: }

    825: typedef struct {
    826:   PetscMPIInt rank;
    827:   DM          dm;
    828:   PetscViewer viewer;
    829:   DMLabel     label;
    830: } LabelCtx;

    832: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
    833: {
    834:   PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
    835:   DMLabel         label  = ((LabelCtx *) op_data)->label;
    836:   IS              stratumIS;
    837:   const PetscInt *ind;
    838:   PetscInt        value, N, i;
    839:   const char     *lname;
    840:   char            group[PETSC_MAX_PATH_LEN];
    841:   PetscErrorCode  ierr;

    843:   PetscOptionsStringToInt(name, &value);
    844:   ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
    845:   PetscObjectSetName((PetscObject) stratumIS, "indices");
    846:   PetscObjectGetName((PetscObject) label, &lname);
    847:   PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
    848:   PetscViewerHDF5PushGroup(viewer, group);
    849:   {
    850:     /* Force serial load */
    851:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
    852:     PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
    853:     PetscLayoutSetSize(stratumIS->map, N);
    854:   }
    855:   ISLoad(stratumIS, viewer);
    856:   PetscViewerHDF5PopGroup(viewer);
    857:   ISGetLocalSize(stratumIS, &N);
    858:   ISGetIndices(stratumIS, &ind);
    859:   for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
    860:   ISRestoreIndices(stratumIS, &ind);
    861:   ISDestroy(&stratumIS);
    862:   return 0;
    863: }

    865: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
    866: {
    867:   DM             dm  = ((LabelCtx *) op_data)->dm;
    868:   hsize_t        idx = 0;
    870:   herr_t         err;

    872:   DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
    873:   DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
    874:   PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
    875:   return err;
    876: }

    878: PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM dm, PetscViewer viewer)
    879: {
    880:   LabelCtx        ctx;
    881:   hid_t           fileId, groupId;
    882:   hsize_t         idx = 0;
    883:   PetscErrorCode  ierr;

    886:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
    887:   ctx.dm     = dm;
    888:   ctx.viewer = viewer;
    889:   PetscViewerHDF5PushGroup(viewer, "/labels");
    890:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
    891:   PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
    892:   PetscStackCallHDF5(H5Gclose,(groupId));
    893:   PetscViewerHDF5PopGroup(viewer);
    894:   return(0);
    895: }

    897: /* The first version will read everything onto proc 0, letting the user distribute
    898:    The next will create a naive partition, and then rebalance after reading
    899: */
    900: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
    901: {
    902:   PetscSection    coordSection;
    903:   Vec             coordinates;
    904:   IS              orderIS, conesIS, cellsIS, orntsIS;
    905:   const PetscInt *order, *cones, *cells, *ornts;
    906:   PetscReal       lengthScale;
    907:   PetscInt       *cone, *ornt;
    908:   PetscInt        dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
    909:   PetscMPIInt     rank;
    910:   PetscErrorCode  ierr;

    913:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
    914:   /* Read toplogy */
    915:   PetscViewerHDF5PushGroup(viewer, "/topology");
    916:   ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
    917:   PetscObjectSetName((PetscObject) orderIS, "order");
    918:   ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
    919:   PetscObjectSetName((PetscObject) conesIS, "cones");
    920:   ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
    921:   PetscObjectSetName((PetscObject) cellsIS, "cells");
    922:   ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
    923:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
    924:   PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
    925:   DMSetDimension(dm, dim);
    926:   {
    927:     /* Force serial load */
    928:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
    929:     PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
    930:     PetscLayoutSetSize(orderIS->map, pEnd);
    931:     PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
    932:     PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
    933:     PetscLayoutSetSize(conesIS->map, pEnd);
    934:     pEnd = !rank ? pEnd : 0;
    935:     PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
    936:     PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
    937:     PetscLayoutSetSize(cellsIS->map, N);
    938:     PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
    939:     PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
    940:     PetscLayoutSetSize(orntsIS->map, N);
    941:   }
    942:   ISLoad(orderIS, viewer);
    943:   ISLoad(conesIS, viewer);
    944:   ISLoad(cellsIS, viewer);
    945:   ISLoad(orntsIS, viewer);
    946:   PetscViewerHDF5PopGroup(viewer);
    947:   /* Read geometry */
    948:   PetscViewerHDF5PushGroup(viewer, "/geometry");
    949:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
    950:   PetscObjectSetName((PetscObject) coordinates, "vertices");
    951:   {
    952:     /* Force serial load */
    953:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
    954:     VecSetSizes(coordinates, !rank ? N : 0, N);
    955:     VecSetBlockSize(coordinates, spatialDim);
    956:   }
    957:   VecLoad(coordinates, viewer);
    958:   PetscViewerHDF5PopGroup(viewer);
    959:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
    960:   VecScale(coordinates, 1.0/lengthScale);
    961:   VecGetLocalSize(coordinates, &numVertices);
    962:   VecGetBlockSize(coordinates, &spatialDim);
    963:   numVertices /= spatialDim;
    964:   /* Create Plex */
    965:   DMPlexSetChart(dm, 0, pEnd);
    966:   ISGetIndices(orderIS, &order);
    967:   ISGetIndices(conesIS, &cones);
    968:   for (p = 0; p < pEnd; ++p) {
    969:     DMPlexSetConeSize(dm, order[p], cones[p]);
    970:     maxConeSize = PetscMax(maxConeSize, cones[p]);
    971:   }
    972:   DMSetUp(dm);
    973:   ISGetIndices(cellsIS, &cells);
    974:   ISGetIndices(orntsIS, &ornts);
    975:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
    976:   for (p = 0, q = 0; p < pEnd; ++p) {
    977:     for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
    978:     DMPlexSetCone(dm, order[p], cone);
    979:     DMPlexSetConeOrientation(dm, order[p], ornt);
    980:   }
    981:   PetscFree2(cone,ornt);
    982:   ISRestoreIndices(orderIS, &order);
    983:   ISRestoreIndices(conesIS, &cones);
    984:   ISRestoreIndices(cellsIS, &cells);
    985:   ISRestoreIndices(orntsIS, &ornts);
    986:   ISDestroy(&orderIS);
    987:   ISDestroy(&conesIS);
    988:   ISDestroy(&cellsIS);
    989:   ISDestroy(&orntsIS);
    990:   DMPlexSymmetrize(dm);
    991:   DMPlexStratify(dm);
    992:   /* Create coordinates */
    993:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    994:   if (numVertices != vEnd - vStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %d does not match number of vertices %d", numVertices, vEnd - vStart);
    995:   DMGetCoordinateSection(dm, &coordSection);
    996:   PetscSectionSetNumFields(coordSection, 1);
    997:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
    998:   PetscSectionSetChart(coordSection, vStart, vEnd);
    999:   for (v = vStart; v < vEnd; ++v) {
    1000:     PetscSectionSetDof(coordSection, v, spatialDim);
    1001:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
    1002:   }
    1003:   PetscSectionSetUp(coordSection);
    1004:   DMSetCoordinates(dm, coordinates);
    1005:   VecDestroy(&coordinates);
    1006:   /* Read Labels */
    1007:   DMPlexLoadLabels_HDF5_Internal(dm, viewer);
    1008:   return(0);
    1009: }
    1010: #endif