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: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);
      8: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer, DMPlexStorageVersion);
      9: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer, const char[], DMPlexStorageVersion);
     10: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);

     12: PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

     14: static PetscErrorCode PetscViewerPrintVersion_Private(PetscViewer viewer, DMPlexStorageVersion version, char str[], size_t len)
     15: {
     16:   PetscFunctionBegin;
     17:   PetscCall(PetscViewerCheckVersion_Private(viewer, version));
     18:   PetscCall(PetscSNPrintf(str, len, "%d.%d.%d", version->major, version->minor, version->subminor));
     19:   PetscFunctionReturn(PETSC_SUCCESS);
     20: }

     22: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer viewer, const char str[], DMPlexStorageVersion *version)
     23: {
     24:   PetscToken           t;
     25:   const char          *ts;
     26:   PetscInt             i;
     27:   PetscInt             ti[3];
     28:   DMPlexStorageVersion v;

     30:   PetscFunctionBegin;
     31:   PetscCall(PetscTokenCreate(str, '.', &t));
     32:   for (i = 0; i < 3; i++) {
     33:     PetscCall(PetscTokenFind(t, &ts));
     34:     PetscCheck(ts, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Malformed version string %s", str);
     35:     PetscCall(PetscOptionsStringToInt(ts, &ti[i]));
     36:   }
     37:   PetscCall(PetscTokenFind(t, &ts));
     38:   PetscCheck(!ts, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Malformed version string %s", str);
     39:   PetscCall(PetscTokenDestroy(&t));
     40:   PetscCall(PetscNew(&v));
     41:   v->major    = (int)ti[0];
     42:   v->minor    = (int)ti[1];
     43:   v->subminor = (int)ti[2];
     44:   PetscCall(PetscViewerCheckVersion_Private(viewer, v));
     45:   *version = v;
     46:   PetscFunctionReturn(PETSC_SUCCESS);
     47: }

     49: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion v)
     50: {
     51:   PetscFunctionBegin;
     52:   PetscCall(PetscObjectContainerCompose((PetscObject)viewer, key, v, PetscCtxDestroyDefault));
     53:   PetscFunctionReturn(PETSC_SUCCESS);
     54: }

     56: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion *v)
     57: {
     58:   PetscContainer cont;

     60:   PetscFunctionBegin;
     61:   PetscCall(PetscObjectQuery((PetscObject)viewer, key, (PetscObject *)&cont));
     62:   *v = NULL;
     63:   if (cont) PetscCall(PetscContainerGetPointer(cont, v));
     64:   PetscFunctionReturn(PETSC_SUCCESS);
     65: }

     67: /*
     68:   Version log:
     69:   1.0.0 legacy version (default if no "dmplex_storage_version" attribute found in file)
     70:   1.1.0 legacy version, but output VIZ by default
     71:   2.0.0 introduce versioning and multiple topologies
     72:   2.1.0 introduce distributions
     73:   3.0.0 new checkpointing format in Firedrake paper
     74:   3.1.0 new format with IS compression
     75: */
     76: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer viewer, DMPlexStorageVersion version)
     77: {
     78:   PetscBool valid = PETSC_FALSE;

     80:   PetscFunctionBegin;
     81:   switch (version->major) {
     82:   case 1:
     83:     switch (version->minor) {
     84:     case 0:
     85:       switch (version->subminor) {
     86:       case 0:
     87:         valid = PETSC_TRUE;
     88:         break;
     89:       }
     90:       break;
     91:     case 1:
     92:       switch (version->subminor) {
     93:       case 0:
     94:         valid = PETSC_TRUE;
     95:         break;
     96:       }
     97:       break;
     98:     }
     99:     break;
    100:   case 2:
    101:     switch (version->minor) {
    102:     case 0:
    103:       switch (version->subminor) {
    104:       case 0:
    105:         valid = PETSC_TRUE;
    106:         break;
    107:       }
    108:       break;
    109:     case 1:
    110:       switch (version->subminor) {
    111:       case 0:
    112:         valid = PETSC_TRUE;
    113:         break;
    114:       }
    115:       break;
    116:     }
    117:     break;
    118:   case 3:
    119:     switch (version->minor) {
    120:     case 0:
    121:       switch (version->subminor) {
    122:       case 0:
    123:         valid = PETSC_TRUE;
    124:         break;
    125:       }
    126:       break;
    127:     case 1:
    128:       switch (version->subminor) {
    129:       case 0:
    130:         valid = PETSC_TRUE;
    131:         break;
    132:       }
    133:       break;
    134:     }
    135:     break;
    136:   }
    137:   PetscCheck(valid, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "DMPlexStorageVersion %d.%d.%d not supported", version->major, version->minor, version->subminor);
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: static inline PetscBool DMPlexStorageVersionEQ(DMPlexStorageVersion version, int major, int minor, int subminor)
    142: {
    143:   return (PetscBool)(version->major == major && version->minor == minor && version->subminor == subminor);
    144: }

    146: static inline PetscBool DMPlexStorageVersionGE(DMPlexStorageVersion version, int major, int minor, int subminor)
    147: {
    148:   return (PetscBool)((version->major == major && version->minor == minor && version->subminor >= subminor) || (version->major == major && version->minor > minor) || version->major > major);
    149: }

    151: /*@C
    152:   PetscViewerHDF5SetDMPlexStorageVersionWriting - Set the storage version for writing

    154:   Logically collective

    156:   Input Parameters:
    157: + viewer  - The `PetscViewer`
    158: - version - The storage format version

    160:   Level: advanced

    162:   Note:
    163:   The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0.

    165: .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5GetDMPlexStorageVersionWriting()`, `PetscViewerHDF5GetDMPlexStorageVersionReading()`, `PetscViewerHDF5SetDMPlexStorageVersionReading()`
    166: @*/
    167: PetscErrorCode PetscViewerHDF5SetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion version)
    168: {
    169:   const char           ATTR_NAME[] = "dmplex_storage_version";
    170:   DMPlexStorageVersion viewerVersion;
    171:   PetscBool            fileHasVersion;
    172:   char                 fileVersion[16], versionStr[16], viewerVersionStr[16];

    174:   PetscFunctionBegin;
    176:   PetscAssertPointer(version, 2);
    177:   PetscCall(PetscViewerPrintVersion_Private(viewer, version, versionStr, 16));
    178:   PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, &viewerVersion));
    179:   if (viewerVersion) {
    180:     PetscBool flg;

    182:     PetscCall(PetscViewerPrintVersion_Private(viewer, viewerVersion, viewerVersionStr, 16));
    183:     PetscCall(PetscStrcmp(versionStr, viewerVersionStr, &flg));
    184:     PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but viewer already has version %s - cannot mix versions", versionStr, viewerVersionStr);
    185:   }

    187:   PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion));
    188:   if (fileHasVersion) {
    189:     PetscBool flg;
    190:     char     *tmp;

    192:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp));
    193:     PetscCall(PetscStrncpy(fileVersion, tmp, sizeof(fileVersion)));
    194:     PetscCall(PetscFree(tmp));
    195:     PetscCall(PetscStrcmp(fileVersion, versionStr, &flg));
    196:     PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", versionStr, fileVersion);
    197:   } else {
    198:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, versionStr));
    199:   }
    200:   PetscCall(PetscNew(&viewerVersion));
    201:   viewerVersion->major    = version->major;
    202:   viewerVersion->minor    = version->minor;
    203:   viewerVersion->subminor = version->subminor;
    204:   PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, viewerVersion));
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: /*@C
    209:   PetscViewerHDF5GetDMPlexStorageVersionWriting - Get the storage version for writing

    211:   Logically collective

    213:   Input Parameter:
    214: . viewer - The `PetscViewer`

    216:   Output Parameter:
    217: . version - The storage format version

    219:   Options Database Keys:
    220: . -dm_plex_view_hdf5_storage_version num - Overrides the storage format version

    222:   Level: advanced

    224:   Note:
    225:   The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0.

    227: .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5SetDMPlexStorageVersionWriting()`, `PetscViewerHDF5GetDMPlexStorageVersionReading()`, `PetscViewerHDF5SetDMPlexStorageVersionReading()`
    228: @*/
    229: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion *version)
    230: {
    231:   const char ATTR_NAME[] = "dmplex_storage_version";
    232:   PetscBool  fileHasVersion;
    233:   char       optVersion[16], fileVersion[16];

    235:   PetscFunctionBegin;
    237:   PetscAssertPointer(version, 2);
    238:   PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, version));
    239:   if (*version) PetscFunctionReturn(PETSC_SUCCESS);

    241:   PetscCall(PetscStrncpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE, sizeof(fileVersion)));
    242:   PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion));
    243:   if (fileHasVersion) {
    244:     char *tmp;

    246:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp));
    247:     PetscCall(PetscStrncpy(fileVersion, tmp, sizeof(fileVersion)));
    248:     PetscCall(PetscFree(tmp));
    249:   }
    250:   PetscCall(PetscStrncpy(optVersion, fileVersion, sizeof(optVersion)));
    251:   PetscObjectOptionsBegin((PetscObject)viewer);
    252:   PetscCall(PetscOptionsString("-dm_plex_view_hdf5_storage_version", "DMPlex HDF5 viewer storage version", NULL, optVersion, optVersion, sizeof(optVersion), NULL));
    253:   PetscOptionsEnd();
    254:   if (!fileHasVersion) {
    255:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, optVersion));
    256:   } else {
    257:     PetscBool flg;

    259:     PetscCall(PetscStrcmp(fileVersion, optVersion, &flg));
    260:     PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", optVersion, fileVersion);
    261:   }
    262:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT));
    263:   PetscCall(PetscViewerParseVersion_Private(viewer, optVersion, version));
    264:   PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, *version));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: /*@C
    269:   PetscViewerHDF5SetDMPlexStorageVersionReading - Set the storage version for reading

    271:   Logically collective

    273:   Input Parameters:
    274: + viewer  - The `PetscViewer`
    275: - version - The storage format version

    277:   Level: advanced

    279:   Note:
    280:   The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0.

    282: .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5GetDMPlexStorageVersionReading()`, `PetscViewerHDF5GetDMPlexStorageVersionWriting()`, `PetscViewerHDF5SetDMPlexStorageVersionWriting()`
    283: @*/
    284: PetscErrorCode PetscViewerHDF5SetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion version)
    285: {
    286:   const char           ATTR_NAME[] = "dmplex_storage_version";
    287:   DMPlexStorageVersion viewerVersion;
    288:   PetscBool            fileHasVersion;
    289:   char                 versionStr[16], viewerVersionStr[16];

    291:   PetscFunctionBegin;
    293:   PetscAssertPointer(version, 2);
    294:   PetscCall(PetscViewerPrintVersion_Private(viewer, version, versionStr, 16));
    295:   PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, &viewerVersion));
    296:   if (viewerVersion) {
    297:     PetscBool flg;

    299:     PetscCall(PetscViewerPrintVersion_Private(viewer, viewerVersion, viewerVersionStr, 16));
    300:     PetscCall(PetscStrcmp(versionStr, viewerVersionStr, &flg));
    301:     PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but viewer already has version %s - cannot mix versions", versionStr, viewerVersionStr);
    302:   }

    304:   PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion));
    305:   if (fileHasVersion) {
    306:     char     *fileVersion;
    307:     PetscBool flg;

    309:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &fileVersion));
    310:     PetscCall(PetscStrcmp(fileVersion, versionStr, &flg));
    311:     PetscCheck(flg, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "User requested DMPlex storage version %s but file already has version %s - cannot mix versions", versionStr, fileVersion);
    312:     PetscCall(PetscFree(fileVersion));
    313:   }
    314:   PetscCall(PetscNew(&viewerVersion));
    315:   viewerVersion->major    = version->major;
    316:   viewerVersion->minor    = version->minor;
    317:   viewerVersion->subminor = version->subminor;
    318:   PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, viewerVersion));
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: /*@C
    323:   PetscViewerHDF5GetDMPlexStorageVersionReading - Get the storage version for reading

    325:   Logically collective

    327:   Input Parameter:
    328: . viewer - The `PetscViewer`

    330:   Output Parameter:
    331: . version - The storage format version

    333:   Options Database Keys:
    334: . -dm_plex_view_hdf5_storage_version num - Overrides the storage format version

    336:   Level: advanced

    338:   Note:
    339:   The version has major, minor, and subminor integers. Parallel operations are only available for version 3.0.0.

    341: .seealso: [](ch_dmbase), `DM`, `PetscViewerHDF5SetDMPlexStorageVersionReading()`, `PetscViewerHDF5GetDMPlexStorageVersionWriting()`, `PetscViewerHDF5SetDMPlexStorageVersionWriting()`
    342: @*/
    343: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion *version)
    344: {
    345:   const char ATTR_NAME[] = "dmplex_storage_version";
    346:   char      *defaultVersion;
    347:   char      *versionString;

    349:   PetscFunctionBegin;
    351:   PetscAssertPointer(version, 2);
    352:   PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, version));
    353:   if (*version) PetscFunctionReturn(PETSC_SUCCESS);

    355:   //TODO string HDF5 attribute handling is terrible and should be redesigned
    356:   PetscCall(PetscStrallocpy(DMPLEX_STORAGE_VERSION_FIRST, &defaultVersion));
    357:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString));
    358:   PetscCall(PetscViewerParseVersion_Private(viewer, versionString, version));
    359:   PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, *version));
    360:   PetscCall(PetscFree(versionString));
    361:   PetscCall(PetscFree(defaultVersion));
    362:   PetscFunctionReturn(PETSC_SUCCESS);
    363: }

    365: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
    366: {
    367:   PetscFunctionBegin;
    368:   if (((PetscObject)dm)->name) PetscCall(PetscObjectGetName((PetscObject)dm, name));
    369:   else *name = "plex";
    370:   PetscFunctionReturn(PETSC_SUCCESS);
    371: }

    373: PetscErrorCode DMSequenceGetLength_HDF5_Internal(DM dm, const char seqname[], PetscInt *seqlen, PetscViewer viewer)
    374: {
    375:   hid_t       file, group, dset, dspace;
    376:   hsize_t     rdim, *dims;
    377:   const char *groupname;
    378:   PetscBool   has;

    380:   PetscFunctionBegin;
    381:   PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &groupname));
    382:   PetscCall(PetscViewerHDF5HasDataset(viewer, seqname, &has));
    383:   PetscCheck(has, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) \"%s\" not stored in group %s", seqname, groupname);

    385:   PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &file, &group));
    386:   PetscCallHDF5Return(dset, H5Dopen2, (group, seqname, H5P_DEFAULT));
    387:   PetscCallHDF5Return(dspace, H5Dget_space, (dset));
    388:   PetscCallHDF5ReturnNoCheck(rdim, H5Sget_simple_extent_dims, (dspace, NULL, NULL));
    389:   PetscCall(PetscMalloc1(rdim, &dims));
    390:   PetscCallHDF5ReturnNoCheck(rdim, H5Sget_simple_extent_dims, (dspace, dims, NULL));
    391:   *seqlen = (PetscInt)dims[0];
    392:   PetscCall(PetscFree(dims));
    393:   PetscCallHDF5(H5Dclose, (dset));
    394:   PetscCallHDF5(H5Gclose, (group));
    395:   PetscFunctionReturn(PETSC_SUCCESS);
    396: }

    398: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char seqname[], PetscInt seqnum, PetscScalar value, PetscViewer viewer)
    399: {
    400:   Vec         stamp;
    401:   PetscMPIInt rank;

    403:   PetscFunctionBegin;
    404:   if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS);
    405:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
    406:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp));
    407:   PetscCall(VecSetBlockSize(stamp, 1));
    408:   PetscCall(PetscObjectSetName((PetscObject)stamp, seqname));
    409:   if (rank == 0) {
    410:     PetscReal timeScale;
    411:     PetscBool istime;

    413:     PetscCall(PetscStrncmp(seqname, "time", 5, &istime));
    414:     if (istime) {
    415:       PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale));
    416:       value *= timeScale;
    417:     }
    418:     PetscCall(VecSetValue(stamp, 0, value, INSERT_VALUES));
    419:   }
    420:   PetscCall(VecAssemblyBegin(stamp));
    421:   PetscCall(VecAssemblyEnd(stamp));
    422:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/"));
    423:   PetscCall(PetscViewerHDF5PushTimestepping(viewer));
    424:   PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */
    425:   PetscCall(VecView(stamp, viewer));
    426:   PetscCall(PetscViewerHDF5PopTimestepping(viewer));
    427:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    428:   PetscCall(VecDestroy(&stamp));
    429:   PetscFunctionReturn(PETSC_SUCCESS);
    430: }

    432: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char seqname[], PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
    433: {
    434:   Vec         stamp;
    435:   PetscMPIInt rank;

    437:   PetscFunctionBegin;
    438:   if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS);
    439:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
    440:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp));
    441:   PetscCall(VecSetBlockSize(stamp, 1));
    442:   PetscCall(PetscObjectSetName((PetscObject)stamp, seqname));
    443:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/"));
    444:   PetscCall(PetscViewerHDF5PushTimestepping(viewer));
    445:   PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */
    446:   PetscCall(VecLoad(stamp, viewer));
    447:   PetscCall(PetscViewerHDF5PopTimestepping(viewer));
    448:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    449:   if (rank == 0) {
    450:     const PetscScalar *a;
    451:     PetscReal          timeScale;
    452:     PetscBool          istime;

    454:     PetscCall(VecGetArrayRead(stamp, &a));
    455:     *value = a[0];
    456:     PetscCall(VecRestoreArrayRead(stamp, &a));
    457:     PetscCall(PetscStrncmp(seqname, "time", 5, &istime));
    458:     if (istime) {
    459:       PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale));
    460:       *value /= timeScale;
    461:     }
    462:   }
    463:   PetscCall(VecDestroy(&stamp));
    464:   PetscFunctionReturn(PETSC_SUCCESS);
    465: }

    467: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
    468: {
    469:   IS              cutcells = NULL;
    470:   const PetscInt *cutc;
    471:   PetscInt        cellHeight, vStart, vEnd, cStart, cEnd, c;

    473:   PetscFunctionBegin;
    474:   if (!cutLabel) PetscFunctionReturn(PETSC_SUCCESS);
    475:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
    476:   PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
    477:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    478:   /* Label vertices that should be duplicated */
    479:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel));
    480:   PetscCall(DMLabelGetStratumIS(cutLabel, 2, &cutcells));
    481:   if (cutcells) {
    482:     PetscInt n;

    484:     PetscCall(ISGetIndices(cutcells, &cutc));
    485:     PetscCall(ISGetLocalSize(cutcells, &n));
    486:     for (c = 0; c < n; ++c) {
    487:       if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
    488:         PetscInt *closure = NULL;
    489:         PetscInt  closureSize, cl, value;

    491:         PetscCall(DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure));
    492:         for (cl = 0; cl < closureSize * 2; cl += 2) {
    493:           if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
    494:             PetscCall(DMLabelGetValue(cutLabel, closure[cl], &value));
    495:             if (value == 1) PetscCall(DMLabelSetValue(*cutVertexLabel, closure[cl], 1));
    496:           }
    497:         }
    498:         PetscCall(DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure));
    499:       }
    500:     }
    501:     PetscCall(ISRestoreIndices(cutcells, &cutc));
    502:   }
    503:   PetscCall(ISDestroy(&cutcells));
    504:   PetscFunctionReturn(PETSC_SUCCESS);
    505: }

    507: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
    508: {
    509:   DMPlexStorageVersion version;
    510:   DM                   dm;
    511:   DM                   dmBC;
    512:   PetscSection         section, sectionGlobal;
    513:   Vec                  gv;
    514:   const char          *name;
    515:   PetscViewerFormat    format;
    516:   PetscInt             seqnum;
    517:   PetscReal            seqval;
    518:   PetscBool            isseq;

    520:   PetscFunctionBegin;
    521:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    522:   PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
    523:   PetscCall(VecGetDM(v, &dm));
    524:   PetscCall(DMGetLocalSection(dm, &section));
    525:   PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
    526:   PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar)seqval, viewer));
    527:   if (seqnum >= 0) {
    528:     PetscCall(PetscViewerHDF5PushTimestepping(viewer));
    529:     PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
    530:   }
    531:   PetscCall(PetscViewerGetFormat(viewer, &format));
    532:   PetscCall(DMGetOutputDM(dm, &dmBC));
    533:   PetscCall(DMGetGlobalSection(dmBC, &sectionGlobal));
    534:   PetscCall(DMGetGlobalVector(dmBC, &gv));
    535:   PetscCall(PetscObjectGetName((PetscObject)v, &name));
    536:   PetscCall(PetscObjectSetName((PetscObject)gv, name));
    537:   PetscCall(DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv));
    538:   PetscCall(DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv));
    539:   PetscCall(PetscObjectTypeCompare((PetscObject)gv, VECSEQ, &isseq));
    540:   if (format == PETSC_VIEWER_HDF5_VIZ) {
    541:     /* Output visualization representation */
    542:     PetscInt numFields, f;
    543:     DMLabel  cutLabel, cutVertexLabel = NULL;

    545:     PetscCall(PetscSectionGetNumFields(section, &numFields));
    546:     PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
    547:     for (f = 0; f < numFields; ++f) {
    548:       Vec                      subv;
    549:       IS                       is;
    550:       const char              *fname, *fgroup, *componentName, *fname_def = "unnamed";
    551:       char                     subname[PETSC_MAX_PATH_LEN];
    552:       PetscInt                 Nc, Nt = 1;
    553:       PetscInt                *pStart, *pEnd;
    554:       PetscViewerVTKFieldType *ft;

    556:       if (DMPlexStorageVersionEQ(version, 1, 1, 0)) PetscCall(DMPlexGetFieldTypes_Internal(dm, section, f, &Nt, &pStart, &pEnd, &ft));
    557:       else {
    558:         PetscCall(PetscMalloc3(Nt, &pStart, Nt, &pEnd, Nt, &ft));
    559:         PetscCall(DMPlexGetFieldType_Internal(dm, section, f, &pStart[0], &pEnd[0], &ft[0]));
    560:       }
    561:       for (PetscInt t = 0; t < Nt; ++t) {
    562:         size_t lname;

    564:         if (ft[t] == PETSC_VTK_INVALID) continue;
    565:         fgroup = (ft[t] == PETSC_VTK_POINT_VECTOR_FIELD) || (ft[t] == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
    566:         PetscCall(PetscSectionGetFieldName(section, f, &fname));
    567:         if (!fname) fname = fname_def;

    569:         if (!t) {
    570:           PetscCall(PetscViewerHDF5PushGroup(viewer, fgroup));
    571:         } else {
    572:           char group[PETSC_MAX_PATH_LEN];

    574:           PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s_%" PetscInt_FMT, fgroup, t));
    575:           PetscCall(PetscViewerHDF5PushGroup(viewer, group));
    576:         }

    578:         if (cutLabel) {
    579:           const PetscScalar *ga;
    580:           PetscScalar       *suba;
    581:           PetscInt           gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0;

    583:           PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
    584:           PetscCall(PetscSectionGetFieldComponents(section, f, &Nc));
    585:           for (PetscInt p = pStart[t]; p < pEnd[t]; ++p) {
    586:             PetscInt gdof, fdof = 0, val;

    588:             PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
    589:             if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
    590:             subSize += fdof;
    591:             PetscCall(DMLabelGetValue(cutVertexLabel, p, &val));
    592:             if (val == 1) extSize += fdof;
    593:           }
    594:           PetscCall(VecCreate(PetscObjectComm((PetscObject)gv), &subv));
    595:           PetscCall(VecSetSizes(subv, subSize + extSize, PETSC_DETERMINE));
    596:           PetscCall(VecSetBlockSize(subv, Nc));
    597:           PetscCall(VecSetType(subv, VECSTANDARD));
    598:           PetscCall(VecGetOwnershipRange(gv, &gstart, NULL));
    599:           PetscCall(VecGetArrayRead(gv, &ga));
    600:           PetscCall(VecGetArray(subv, &suba));
    601:           for (PetscInt p = pStart[t]; p < pEnd[t]; ++p) {
    602:             PetscInt gdof, goff, val;

    604:             PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
    605:             if (gdof > 0) {
    606:               PetscInt fdof, fc, f2, poff = 0;

    608:               PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
    609:               /* Can get rid of this loop by storing field information in the global section */
    610:               for (f2 = 0; f2 < f; ++f2) {
    611:                 PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
    612:                 poff += fdof;
    613:               }
    614:               PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
    615:               for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff + poff + fc - gstart];
    616:               PetscCall(DMLabelGetValue(cutVertexLabel, p, &val));
    617:               if (val == 1) {
    618:                 for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize + newOff] = ga[goff + poff + fc - gstart];
    619:               }
    620:             }
    621:           }
    622:           PetscCall(VecRestoreArrayRead(gv, &ga));
    623:           PetscCall(VecRestoreArray(subv, &suba));
    624:           PetscCall(DMLabelDestroy(&cutVertexLabel));
    625:         } else {
    626:           PetscCall(PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart[t], pEnd[t], &is, &subv));
    627:         }
    628:         PetscCall(PetscStrlen(name, &lname));
    629:         if (lname) {
    630:           PetscCall(PetscStrncpy(subname, name, sizeof(subname)));
    631:           PetscCall(PetscStrlcat(subname, "_", sizeof(subname)));
    632:         }
    633:         PetscCall(PetscStrlcat(subname, fname, sizeof(subname)));
    634:         PetscCall(PetscObjectSetName((PetscObject)subv, subname));
    635:         if (isseq) PetscCall(VecView_Seq(subv, viewer));
    636:         else PetscCall(VecView_MPI(subv, viewer));
    637:         if (ft[t] == PETSC_VTK_POINT_VECTOR_FIELD || ft[t] == PETSC_VTK_CELL_VECTOR_FIELD) {
    638:           PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "vector"));
    639:         } else {
    640:           PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "scalar"));
    641:         }

    643:         /* Output the component names in the field if available */
    644:         PetscCall(PetscSectionGetFieldComponents(section, f, &Nc));
    645:         for (PetscInt c = 0; c < Nc; ++c) {
    646:           char componentNameLabel[PETSC_MAX_PATH_LEN];
    647:           PetscCall(PetscSectionGetComponentName(section, f, c, &componentName));
    648:           PetscCall(PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%" PetscInt_FMT, c));
    649:           PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, componentNameLabel, PETSC_STRING, componentName));
    650:         }

    652:         if (cutLabel) PetscCall(VecDestroy(&subv));
    653:         else PetscCall(PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart[t], pEnd[t], &is, &subv));
    654:         PetscCall(PetscViewerHDF5PopGroup(viewer));
    655:       }
    656:       if (!DMPlexStorageVersionEQ(version, 1, 1, 0)) PetscCall(PetscFree3(pStart, pEnd, ft));
    657:     }
    658:   } else {
    659:     /* Output full vector */
    660:     PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
    661:     if (isseq) PetscCall(VecView_Seq(gv, viewer));
    662:     else PetscCall(VecView_MPI(gv, viewer));
    663:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    664:   }
    665:   if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
    666:   PetscCall(DMRestoreGlobalVector(dmBC, &gv));
    667:   PetscFunctionReturn(PETSC_SUCCESS);
    668: }

    670: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
    671: {
    672:   DMPlexStorageVersion version;
    673:   DM                   dm;
    674:   Vec                  locv;
    675:   PetscObject          isZero;
    676:   const char          *name;
    677:   PetscReal            time;

    679:   PetscFunctionBegin;
    680:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    681:   PetscCall(PetscInfo(v, "Writing Vec %s storage version %d.%d.%d\n", v->hdr.name, version->major, version->minor, version->subminor));

    683:   PetscCall(VecGetDM(v, &dm));
    684:   PetscCall(DMGetLocalVector(dm, &locv));
    685:   PetscCall(PetscObjectGetName((PetscObject)v, &name));
    686:   PetscCall(PetscObjectSetName((PetscObject)locv, name));
    687:   PetscCall(PetscObjectQuery((PetscObject)v, "__Vec_bc_zero__", &isZero));
    688:   PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", isZero));
    689:   PetscCall(DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv));
    690:   PetscCall(DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv));
    691:   PetscCall(DMGetOutputSequenceNumber(dm, NULL, &time));
    692:   PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL));
    693:   PetscCall(VecView_Plex_Local_HDF5_Internal(locv, viewer));
    694:   if (DMPlexStorageVersionEQ(version, 1, 1, 0)) {
    695:     PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
    696:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
    697:     PetscCall(VecView_Plex_Local_HDF5_Internal(locv, viewer));
    698:     PetscCall(PetscViewerPopFormat(viewer));
    699:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    700:   }
    701:   PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", NULL));
    702:   PetscCall(DMRestoreLocalVector(dm, &locv));
    703:   PetscFunctionReturn(PETSC_SUCCESS);
    704: }

    706: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
    707: {
    708:   PetscBool isseq;

    710:   PetscFunctionBegin;
    711:   PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
    712:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
    713:   if (isseq) PetscCall(VecView_Seq(v, viewer));
    714:   else PetscCall(VecView_MPI(v, viewer));
    715:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    716:   PetscFunctionReturn(PETSC_SUCCESS);
    717: }

    719: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
    720: {
    721:   DM          dm;
    722:   Vec         locv;
    723:   const char *name;
    724:   PetscInt    seqnum;

    726:   PetscFunctionBegin;
    727:   PetscCall(VecGetDM(v, &dm));
    728:   PetscCall(DMGetLocalVector(dm, &locv));
    729:   PetscCall(PetscObjectGetName((PetscObject)v, &name));
    730:   PetscCall(PetscObjectSetName((PetscObject)locv, name));
    731:   PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL));
    732:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
    733:   if (seqnum >= 0) {
    734:     PetscCall(PetscViewerHDF5PushTimestepping(viewer));
    735:     PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
    736:   }
    737:   PetscCall(VecLoad_Plex_Local(locv, viewer));
    738:   if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
    739:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    740:   PetscCall(DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v));
    741:   PetscCall(DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v));
    742:   PetscCall(DMRestoreLocalVector(dm, &locv));
    743:   PetscFunctionReturn(PETSC_SUCCESS);
    744: }

    746: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
    747: {
    748:   DM       dm;
    749:   PetscInt seqnum;

    751:   PetscFunctionBegin;
    752:   PetscCall(VecGetDM(v, &dm));
    753:   PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL));
    754:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
    755:   if (seqnum >= 0) {
    756:     PetscCall(PetscViewerHDF5PushTimestepping(viewer));
    757:     PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
    758:   }
    759:   PetscCall(VecLoad_Default(v, viewer));
    760:   if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
    761:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    762:   PetscFunctionReturn(PETSC_SUCCESS);
    763: }

    765: static PetscErrorCode DMPlexDistributionView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
    766: {
    767:   DMPlexStorageVersion version;
    768:   MPI_Comm             comm;
    769:   PetscMPIInt          size, rank;
    770:   PetscInt             size_petsc_int;
    771:   const char          *topologydm_name, *distribution_name;
    772:   const PetscInt      *gpoint;
    773:   PetscInt             pStart, pEnd, p;
    774:   PetscSF              pointSF;
    775:   PetscInt             nroots, nleaves;
    776:   const PetscInt      *ilocal;
    777:   const PetscSFNode   *iremote;
    778:   IS                   chartSizesIS, ownersIS, gpointsIS;
    779:   PetscInt            *chartSize, *owners, *gpoints;
    780:   PetscBool            ocompress;

    782:   PetscFunctionBegin;
    783:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    784:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    785:   PetscCallMPI(MPI_Comm_size(comm, &size));
    786:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    787:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    788:   PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
    789:   if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
    790:   PetscCall(PetscLogEventBegin(DMPLEX_DistributionView, viewer, 0, 0, 0));
    791:   size_petsc_int = (PetscInt)size;
    792:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size_petsc_int));
    793:   PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
    794:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    795:   PetscCall(PetscMalloc1(1, &chartSize));
    796:   *chartSize = pEnd - pStart;
    797:   PetscCall(PetscMalloc1(*chartSize, &owners));
    798:   PetscCall(PetscMalloc1(*chartSize, &gpoints));
    799:   PetscCall(DMGetPointSF(dm, &pointSF));
    800:   PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
    801:   for (p = pStart; p < pEnd; ++p) {
    802:     PetscInt gp = gpoint[p - pStart];

    804:     owners[p - pStart]  = rank;
    805:     gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp);
    806:   }
    807:   for (p = 0; p < nleaves; ++p) {
    808:     PetscInt ilocalp = (ilocal ? ilocal[p] : p);

    810:     owners[ilocalp] = iremote[p].rank;
    811:   }
    812:   PetscCall(ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS));
    813:   PetscCall(ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS));
    814:   PetscCall(ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS));
    815:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
    816:     PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
    817:     PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
    818:   }
    819:   PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
    820:   PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
    821:   PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
    822:   PetscCall(ISView(chartSizesIS, viewer));
    823:   PetscCall(ISView(ownersIS, viewer));
    824:   PetscCall(ISView(gpointsIS, viewer));
    825:   PetscCall(ISDestroy(&chartSizesIS));
    826:   PetscCall(ISDestroy(&ownersIS));
    827:   PetscCall(ISDestroy(&gpointsIS));
    828:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
    829:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
    830:   PetscCall(PetscLogEventEnd(DMPLEX_DistributionView, viewer, 0, 0, 0));
    831:   PetscFunctionReturn(PETSC_SUCCESS);
    832: }

    834: static PetscErrorCode DMPlexTopologyView_HDF5_Inner_Private(DM dm, IS globalPointNumbers, PetscViewer viewer, PetscInt pStart, PetscInt pEnd, const char pointsName[], const char coneSizesName[], const char conesName[], const char orientationsName[])
    835: {
    836:   DMPlexStorageVersion version;
    837:   IS                   coneSizesIS, conesIS, orientationsIS;
    838:   PetscInt            *coneSizes, *cones, *orientations;
    839:   const PetscInt      *gpoint;
    840:   PetscInt             nPoints = 0, conesSize = 0;
    841:   PetscInt             p, c, s;
    842:   PetscBool            ocompress;
    843:   MPI_Comm             comm;

    845:   PetscFunctionBegin;
    846:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    847:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    848:   PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
    849:   for (p = pStart; p < pEnd; ++p) {
    850:     if (gpoint[p] >= 0) {
    851:       PetscInt coneSize;

    853:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    854:       nPoints += 1;
    855:       conesSize += coneSize;
    856:     }
    857:   }
    858:   PetscCall(PetscMalloc1(nPoints, &coneSizes));
    859:   PetscCall(PetscMalloc1(conesSize, &cones));
    860:   PetscCall(PetscMalloc1(conesSize, &orientations));
    861:   for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
    862:     if (gpoint[p] >= 0) {
    863:       const PetscInt *cone, *ornt;
    864:       PetscInt        coneSize, cp;

    866:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
    867:       PetscCall(DMPlexGetCone(dm, p, &cone));
    868:       PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
    869:       coneSizes[s] = coneSize;
    870:       for (cp = 0; cp < coneSize; ++cp, ++c) {
    871:         cones[c]        = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]];
    872:         orientations[c] = ornt[cp];
    873:       }
    874:       ++s;
    875:     }
    876:   }
    877:   PetscCheck(s == nPoints, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %" PetscInt_FMT " != %" PetscInt_FMT, s, nPoints);
    878:   PetscCheck(c == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %" PetscInt_FMT " != %" PetscInt_FMT, c, conesSize);
    879:   PetscCall(ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS));
    880:   PetscCall(ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS));
    881:   PetscCall(ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS));
    882:   PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
    883:   PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
    884:   PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
    885:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
    886:     PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
    887:     PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
    888:   }
    889:   PetscCall(ISView(coneSizesIS, viewer));
    890:   PetscCall(ISView(conesIS, viewer));
    891:   PetscCall(ISView(orientationsIS, viewer));
    892:   PetscCall(ISDestroy(&coneSizesIS));
    893:   PetscCall(ISDestroy(&conesIS));
    894:   PetscCall(ISDestroy(&orientationsIS));
    895:   if (pointsName) {
    896:     IS        pointsIS;
    897:     PetscInt *points;

    899:     PetscCall(PetscMalloc1(nPoints, &points));
    900:     for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
    901:       if (gpoint[p] >= 0) {
    902:         points[s] = gpoint[p];
    903:         ++s;
    904:       }
    905:     }
    906:     PetscCall(ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS));
    907:     PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
    908:     PetscCall(ISView(pointsIS, viewer));
    909:     PetscCall(ISDestroy(&pointsIS));
    910:   }
    911:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
    912:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
    913:   PetscFunctionReturn(PETSC_SUCCESS);
    914: }

    916: static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
    917: {
    918:   const char *pointsName, *coneSizesName, *conesName, *orientationsName;
    919:   PetscInt    pStart, pEnd, dim;

    921:   PetscFunctionBegin;
    922:   pointsName       = "order";
    923:   coneSizesName    = "cones";
    924:   conesName        = "cells";
    925:   orientationsName = "orientation";
    926:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    927:   PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName));
    928:   PetscCall(DMGetDimension(dm, &dim));
    929:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim));
    930:   PetscFunctionReturn(PETSC_SUCCESS);
    931: }

    933: //TODO get this numbering right away without needing this function
    934: /* Renumber global point numbers so that they are 0-based per stratum */
    935: static PetscErrorCode RenumberGlobalPointNumbersPerStratum_Private(DM dm, IS globalPointNumbers, IS *newGlobalPointNumbers, IS *strataPermutation)
    936: {
    937:   PetscInt        d, depth, p, n;
    938:   PetscInt       *offsets;
    939:   const PetscInt *gpn;
    940:   PetscInt       *ngpn;
    941:   MPI_Comm        comm;

    943:   PetscFunctionBegin;
    944:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    945:   PetscCall(ISGetLocalSize(globalPointNumbers, &n));
    946:   PetscCall(ISGetIndices(globalPointNumbers, &gpn));
    947:   PetscCall(PetscMalloc1(n, &ngpn));
    948:   PetscCall(DMPlexGetDepth(dm, &depth));
    949:   PetscCall(PetscMalloc1(depth + 1, &offsets));
    950:   for (d = 0; d <= depth; d++) {
    951:     PetscInt pStart, pEnd;

    953:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    954:     offsets[d] = PETSC_INT_MAX;
    955:     for (p = pStart; p < pEnd; p++) {
    956:       if (gpn[p] >= 0 && gpn[p] < offsets[d]) offsets[d] = gpn[p];
    957:     }
    958:   }
    959:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, offsets, depth + 1, MPIU_INT, MPI_MIN, comm));
    960:   for (d = 0; d <= depth; d++) {
    961:     PetscInt pStart, pEnd;

    963:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    964:     for (p = pStart; p < pEnd; p++) ngpn[p] = gpn[p] - PetscSign(gpn[p]) * offsets[d];
    965:   }
    966:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpn));
    967:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)globalPointNumbers), n, ngpn, PETSC_OWN_POINTER, newGlobalPointNumbers));
    968:   {
    969:     PetscInt *perm;

    971:     PetscCall(PetscMalloc1(depth + 1, &perm));
    972:     for (d = 0; d <= depth; d++) perm[d] = d;
    973:     PetscCall(PetscSortIntWithPermutation(depth + 1, offsets, perm));
    974:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, depth + 1, perm, PETSC_OWN_POINTER, strataPermutation));
    975:   }
    976:   PetscCall(PetscFree(offsets));
    977:   PetscFunctionReturn(PETSC_SUCCESS);
    978: }

    980: static PetscErrorCode DMPlexTopologyView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
    981: {
    982:   IS          globalPointNumbers0, strataPermutation;
    983:   const char *coneSizesName, *conesName, *orientationsName;
    984:   PetscInt    depth, d;
    985:   MPI_Comm    comm;

    987:   PetscFunctionBegin;
    988:   coneSizesName    = "cone_sizes";
    989:   conesName        = "cones";
    990:   orientationsName = "orientations";
    991:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    992:   PetscCall(DMPlexGetDepth(dm, &depth));
    993:   {
    994:     PetscInt dim;

    996:     PetscCall(DMGetDimension(dm, &dim));
    997:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "cell_dim", PETSC_INT, &dim));
    998:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "depth", PETSC_INT, &depth));
    999:   }

    1001:   PetscCall(RenumberGlobalPointNumbersPerStratum_Private(dm, globalPointNumbers, &globalPointNumbers0, &strataPermutation));
    1002:   /* TODO dirty trick to save serial IS using the same parallel viewer */
    1003:   {
    1004:     IS              spOnComm;
    1005:     PetscInt        n   = 0, N;
    1006:     const PetscInt *idx = NULL;
    1007:     const PetscInt *old;
    1008:     PetscMPIInt     rank;

    1010:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1011:     PetscCall(ISGetLocalSize(strataPermutation, &N));
    1012:     PetscCall(ISGetIndices(strataPermutation, &old));
    1013:     if (!rank) {
    1014:       n   = N;
    1015:       idx = old;
    1016:     }
    1017:     PetscCall(ISCreateGeneral(comm, n, idx, PETSC_COPY_VALUES, &spOnComm));
    1018:     PetscCall(ISRestoreIndices(strataPermutation, &old));
    1019:     PetscCall(ISDestroy(&strataPermutation));
    1020:     strataPermutation = spOnComm;
    1021:   }
    1022:   PetscCall(PetscObjectSetName((PetscObject)strataPermutation, "permutation"));
    1023:   PetscCall(ISView(strataPermutation, viewer));
    1024:   PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
    1025:   for (d = 0; d <= depth; d++) {
    1026:     PetscInt pStart, pEnd;
    1027:     char     group[128];

    1029:     PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, d));
    1030:     PetscCall(PetscViewerHDF5PushGroup(viewer, group));
    1031:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
    1032:     PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers0, viewer, pStart, pEnd, NULL, coneSizesName, conesName, orientationsName));
    1033:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    1034:   }
    1035:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */
    1036:   PetscCall(ISDestroy(&globalPointNumbers0));
    1037:   PetscCall(ISDestroy(&strataPermutation));
    1038:   PetscFunctionReturn(PETSC_SUCCESS);
    1039: }

    1041: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
    1042: {
    1043:   DMPlexStorageVersion version;
    1044:   const char          *topologydm_name;
    1045:   char                 group[PETSC_MAX_PATH_LEN];

    1047:   PetscFunctionBegin;
    1048:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    1049:   PetscCall(PetscInfo(dm, "Writing DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor));
    1050:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1051:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
    1052:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
    1053:   } else {
    1054:     PetscCall(PetscStrncpy(group, "/", sizeof(group)));
    1055:   }
    1056:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));

    1058:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
    1059:   if (version->major < 3) {
    1060:     PetscCall(DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer));
    1061:   } else {
    1062:     /* since DMPlexStorageVersion 3.0.0 */
    1063:     PetscCall(DMPlexTopologyView_HDF5_Private(dm, globalPointNumbers, viewer));
    1064:   }
    1065:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */

    1067:   if (DMPlexStorageVersionGE(version, 2, 1, 0)) {
    1068:     const char *distribution_name;

    1070:     PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
    1071:     PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
    1072:     PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
    1073:     PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
    1074:     PetscCall(DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer));
    1075:     PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
    1076:     PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
    1077:   }

    1079:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1080:   PetscFunctionReturn(PETSC_SUCCESS);
    1081: }

    1083: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
    1084: {
    1085:   PetscSF         sfPoint;
    1086:   DMLabel         cutLabel, cutVertexLabel         = NULL;
    1087:   IS              globalVertexNumbers, cutvertices = NULL;
    1088:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
    1089:   PetscInt       *vertices;
    1090:   PetscInt        conesSize = 0;
    1091:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;

    1093:   PetscFunctionBegin;
    1094:   *numCorners = 0;
    1095:   PetscCall(DMGetDimension(dm, &dim));
    1096:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    1097:   PetscCall(ISGetIndices(globalCellNumbers, &gcell));

    1099:   for (cell = cStart; cell < cEnd; ++cell) {
    1100:     PetscInt *closure = NULL;
    1101:     PetscInt  closureSize, v, Nc = 0;

    1103:     if (gcell[cell] < 0) continue;
    1104:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
    1105:     for (v = 0; v < closureSize * 2; v += 2) {
    1106:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
    1107:     }
    1108:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
    1109:     conesSize += Nc;
    1110:     if (!numCornersLocal) numCornersLocal = Nc;
    1111:     else if (numCornersLocal != Nc) numCornersLocal = 1;
    1112:   }
    1113:   PetscCallMPI(MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
    1114:   PetscCheck(!numCornersLocal || !(numCornersLocal != *numCorners || *numCorners == 0), PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
    1115:   /* Handle periodic cuts by identifying vertices which should be duplicated */
    1116:   PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
    1117:   PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
    1118:   if (cutVertexLabel) PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices));
    1119:   if (cutvertices) {
    1120:     PetscCall(ISGetIndices(cutvertices, &cutverts));
    1121:     PetscCall(ISGetLocalSize(cutvertices, &vExtra));
    1122:   }
    1123:   PetscCall(DMGetPointSF(dm, &sfPoint));
    1124:   if (cutLabel) {
    1125:     const PetscInt    *ilocal;
    1126:     const PetscSFNode *iremote;
    1127:     PetscInt           nroots, nleaves;

    1129:     PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote));
    1130:     if (nleaves < 0) {
    1131:       PetscCall(PetscObjectReference((PetscObject)sfPoint));
    1132:     } else {
    1133:       PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint));
    1134:       PetscCall(PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES));
    1135:     }
    1136:   } else {
    1137:     PetscCall(PetscObjectReference((PetscObject)sfPoint));
    1138:   }
    1139:   /* Number all vertices */
    1140:   PetscCall(DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers));
    1141:   PetscCall(PetscSFDestroy(&sfPoint));
    1142:   /* Create cones */
    1143:   PetscCall(ISGetIndices(globalVertexNumbers, &gvertex));
    1144:   PetscCall(PetscMalloc1(conesSize, &vertices));
    1145:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
    1146:     PetscInt *closure = NULL;
    1147:     PetscInt  closureSize, Nc = 0, p, value = -1;
    1148:     PetscBool replace;

    1150:     if (gcell[cell] < 0) continue;
    1151:     if (cutLabel) PetscCall(DMLabelGetValue(cutLabel, cell, &value));
    1152:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
    1153:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
    1154:     for (p = 0; p < closureSize * 2; p += 2) {
    1155:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p];
    1156:     }
    1157:     PetscCall(DMPlexReorderCell(dm, cell, closure));
    1158:     for (p = 0; p < Nc; ++p) {
    1159:       PetscInt nv, gv = gvertex[closure[p] - vStart];

    1161:       if (replace) {
    1162:         PetscCall(PetscFindInt(closure[p], vExtra, cutverts, &nv));
    1163:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
    1164:       }
    1165:       vertices[v++] = gv < 0 ? -(gv + 1) : gv;
    1166:     }
    1167:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
    1168:   }
    1169:   PetscCall(ISRestoreIndices(globalVertexNumbers, &gvertex));
    1170:   PetscCall(ISDestroy(&globalVertexNumbers));
    1171:   PetscCall(ISRestoreIndices(globalCellNumbers, &gcell));
    1172:   if (cutvertices) PetscCall(ISRestoreIndices(cutvertices, &cutverts));
    1173:   PetscCall(ISDestroy(&cutvertices));
    1174:   PetscCall(DMLabelDestroy(&cutVertexLabel));
    1175:   PetscCheck(v == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %" PetscInt_FMT " != %" PetscInt_FMT, v, conesSize);
    1176:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS));
    1177:   PetscCall(PetscLayoutSetBlockSize((*cellIS)->map, *numCorners));
    1178:   PetscCall(PetscObjectSetName((PetscObject)*cellIS, "cells"));
    1179:   PetscFunctionReturn(PETSC_SUCCESS);
    1180: }

    1182: static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer)
    1183: {
    1184:   DM       cdm;
    1185:   DMLabel  depthLabel, ctLabel;
    1186:   IS       cellIS;
    1187:   PetscInt dim, depth, cellHeight, c, n = 0;

    1189:   PetscFunctionBegin;
    1190:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
    1191:   PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));

    1193:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1194:   PetscCall(DMGetDimension(dm, &dim));
    1195:   PetscCall(DMPlexGetDepth(dm, &depth));
    1196:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1197:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
    1198:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
    1199:   PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
    1200:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
    1201:     const DMPolytopeType ict = (DMPolytopeType)c;
    1202:     PetscInt             pStart, pEnd, dep, numCorners;
    1203:     PetscBool            output = PETSC_FALSE, doOutput;

    1205:     if (ict == DM_POLYTOPE_FV_GHOST) continue;
    1206:     PetscCall(DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd));
    1207:     if (pStart >= 0) {
    1208:       PetscCall(DMLabelGetValue(depthLabel, pStart, &dep));
    1209:       if (dep == depth - cellHeight) output = PETSC_TRUE;
    1210:     }
    1211:     PetscCallMPI(MPIU_Allreduce(&output, &doOutput, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
    1212:     if (!doOutput) continue;
    1213:     PetscCall(CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS));
    1214:     if (!n) {
    1215:       PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/topology"));
    1216:     } else {
    1217:       char group[PETSC_MAX_PATH_LEN];

    1219:       PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n));
    1220:       PetscCall(PetscViewerHDF5PushGroup(viewer, group));
    1221:     }
    1222:     PetscCall(ISView(cellIS, viewer));
    1223:     PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners));
    1224:     PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim));
    1225:     PetscCall(ISDestroy(&cellIS));
    1226:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    1227:     ++n;
    1228:   }
    1229:   PetscFunctionReturn(PETSC_SUCCESS);
    1230: }

    1232: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
    1233: {
    1234:   DM        cdm;
    1235:   Vec       coordinates, newcoords;
    1236:   PetscReal lengthScale;
    1237:   PetscInt  m, M, bs;

    1239:   PetscFunctionBegin;
    1240:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    1241:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1242:   PetscCall(DMGetCoordinates(dm, &coordinates));
    1243:   PetscCall(VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords));
    1244:   PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
    1245:   PetscCall(VecGetSize(coordinates, &M));
    1246:   PetscCall(VecGetLocalSize(coordinates, &m));
    1247:   PetscCall(VecSetSizes(newcoords, m, M));
    1248:   PetscCall(VecGetBlockSize(coordinates, &bs));
    1249:   PetscCall(VecSetBlockSize(newcoords, bs));
    1250:   PetscCall(VecSetType(newcoords, VECSTANDARD));
    1251:   PetscCall(VecCopy(coordinates, newcoords));
    1252:   PetscCall(VecScale(newcoords, lengthScale));
    1253:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
    1254:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
    1255:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
    1256:   PetscCall(VecView(newcoords, viewer));
    1257:   PetscCall(PetscViewerPopFormat(viewer));
    1258:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1259:   PetscCall(VecDestroy(&newcoords));
    1260:   PetscFunctionReturn(PETSC_SUCCESS);
    1261: }

    1263: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
    1264: {
    1265:   DM                   cdm;
    1266:   Vec                  coords, newcoords;
    1267:   PetscInt             m, M, bs;
    1268:   PetscReal            lengthScale;
    1269:   PetscBool            viewSection = PETSC_TRUE, ocompress;
    1270:   const char          *topologydm_name, *coordinatedm_name, *coordinates_name;
    1271:   PetscViewerFormat    format;
    1272:   DMPlexStorageVersion version;

    1274:   PetscFunctionBegin;
    1275:   PetscCall(PetscViewerGetFormat(viewer, &format));
    1276:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    1277:   if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) {
    1278:     PetscCall(DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer));
    1279:     PetscFunctionReturn(PETSC_SUCCESS);
    1280:   }
    1281:   /* since 2.0.0 */
    1282:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
    1283:     PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
    1284:     PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
    1285:   }
    1286:   PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_view_coordinate_section", &viewSection, NULL));
    1287:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1288:   PetscCall(DMGetCoordinates(dm, &coords));
    1289:   PetscCall(PetscObjectGetName((PetscObject)cdm, &coordinatedm_name));
    1290:   PetscCall(PetscObjectGetName((PetscObject)coords, &coordinates_name));
    1291:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1292:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    1293:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    1294:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name));
    1295:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name));
    1296:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1297:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1298:   if (viewSection) PetscCall(DMPlexSectionView(dm, viewer, cdm));
    1299:   PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &newcoords));
    1300:   PetscCall(PetscObjectSetName((PetscObject)newcoords, coordinates_name));
    1301:   PetscCall(VecGetSize(coords, &M));
    1302:   PetscCall(VecGetLocalSize(coords, &m));
    1303:   PetscCall(VecSetSizes(newcoords, m, M));
    1304:   PetscCall(VecGetBlockSize(coords, &bs));
    1305:   PetscCall(VecSetBlockSize(newcoords, bs));
    1306:   PetscCall(VecSetType(newcoords, VECSTANDARD));
    1307:   PetscCall(VecCopy(coords, newcoords));
    1308:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    1309:   PetscCall(VecScale(newcoords, lengthScale));
    1310:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
    1311:   PetscCall(DMPlexGlobalVectorView(dm, viewer, cdm, newcoords));
    1312:   PetscCall(PetscViewerPopFormat(viewer));
    1313:   PetscCall(VecDestroy(&newcoords));
    1314:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
    1315:   PetscFunctionReturn(PETSC_SUCCESS);
    1316: }

    1318: static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer)
    1319: {
    1320:   DM               cdm;
    1321:   Vec              coordinatesLocal, newcoords;
    1322:   PetscSection     cSection, cGlobalSection;
    1323:   PetscScalar     *coords, *ncoords;
    1324:   DMLabel          cutLabel, cutVertexLabel = NULL;
    1325:   const PetscReal *L;
    1326:   PetscReal        lengthScale;
    1327:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
    1328:   PetscBool        localized, embedded;

    1330:   PetscFunctionBegin;
    1331:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    1332:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    1333:   PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
    1334:   PetscCall(VecGetBlockSize(coordinatesLocal, &bs));
    1335:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
    1336:   if (localized == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
    1337:   PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
    1338:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1339:   PetscCall(DMGetLocalSection(cdm, &cSection));
    1340:   PetscCall(DMGetGlobalSection(cdm, &cGlobalSection));
    1341:   PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
    1342:   N = 0;

    1344:   PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
    1345:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &newcoords));
    1346:   PetscCall(PetscSectionGetDof(cSection, vStart, &dof));
    1347:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof));
    1348:   embedded = (PetscBool)(L && dof == 2 && !cutLabel);
    1349:   if (cutVertexLabel) {
    1350:     PetscCall(DMLabelGetStratumSize(cutVertexLabel, 1, &v));
    1351:     N += dof * v;
    1352:   }
    1353:   for (v = vStart; v < vEnd; ++v) {
    1354:     PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
    1355:     if (dof < 0) continue;
    1356:     if (embedded) N += dof + 1;
    1357:     else N += dof;
    1358:   }
    1359:   if (embedded) PetscCall(VecSetBlockSize(newcoords, bs + 1));
    1360:   else PetscCall(VecSetBlockSize(newcoords, bs));
    1361:   PetscCall(VecSetSizes(newcoords, N, PETSC_DETERMINE));
    1362:   PetscCall(VecSetType(newcoords, VECSTANDARD));
    1363:   PetscCall(VecGetArray(coordinatesLocal, &coords));
    1364:   PetscCall(VecGetArray(newcoords, &ncoords));
    1365:   coordSize = 0;
    1366:   for (v = vStart; v < vEnd; ++v) {
    1367:     PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
    1368:     PetscCall(PetscSectionGetOffset(cSection, v, &off));
    1369:     if (dof < 0) continue;
    1370:     if (embedded) {
    1371:       if (L && (L[0] > 0.0) && (L[1] > 0.0)) {
    1372:         PetscReal theta, phi, r, R;
    1373:         /* XY-periodic */
    1374:         /* Suppose its an y-z circle, then
    1375:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
    1376:            and the circle in that plane is
    1377:              \hat r cos(phi) + \hat x sin(phi) */
    1378:         theta                = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1];
    1379:         phi                  = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0];
    1380:         r                    = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]);
    1381:         R                    = L[1] / (2.0 * PETSC_PI);
    1382:         ncoords[coordSize++] = PetscSinReal(phi) * r;
    1383:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
    1384:         ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
    1385:       } else if (L && (L[0] > 0.0)) {
    1386:         /* X-periodic */
    1387:         ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
    1388:         ncoords[coordSize++] = coords[off + 1];
    1389:         ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
    1390:       } else if (L && (L[1] > 0.0)) {
    1391:         /* Y-periodic */
    1392:         ncoords[coordSize++] = coords[off + 0];
    1393:         ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
    1394:         ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
    1395: #if 0
    1396:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
    1397:         PetscReal phi, r, R;
    1398:         /* Mobius strip */
    1399:         /* Suppose its an x-z circle, then
    1400:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
    1401:            and in that plane we rotate by pi as we go around the circle
    1402:              \hat r cos(phi/2) + \hat y sin(phi/2) */
    1403:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
    1404:         R     = L[0];
    1405:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
    1406:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
    1407:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
    1408:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
    1409: #endif
    1410:       } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
    1411:     } else {
    1412:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d];
    1413:     }
    1414:   }
    1415:   if (cutVertexLabel) {
    1416:     IS              vertices;
    1417:     const PetscInt *verts;
    1418:     PetscInt        n;

    1420:     PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &vertices));
    1421:     if (vertices) {
    1422:       PetscCall(ISGetIndices(vertices, &verts));
    1423:       PetscCall(ISGetLocalSize(vertices, &n));
    1424:       for (v = 0; v < n; ++v) {
    1425:         PetscCall(PetscSectionGetDof(cSection, verts[v], &dof));
    1426:         PetscCall(PetscSectionGetOffset(cSection, verts[v], &off));
    1427:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0);
    1428:       }
    1429:       PetscCall(ISRestoreIndices(vertices, &verts));
    1430:       PetscCall(ISDestroy(&vertices));
    1431:     }
    1432:   }
    1433:   PetscCheck(coordSize == N, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %" PetscInt_FMT " != %" PetscInt_FMT, coordSize, N);
    1434:   PetscCall(DMLabelDestroy(&cutVertexLabel));
    1435:   PetscCall(VecRestoreArray(coordinatesLocal, &coords));
    1436:   PetscCall(VecRestoreArray(newcoords, &ncoords));
    1437:   PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
    1438:   PetscCall(VecScale(newcoords, lengthScale));
    1439:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
    1440:   PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
    1441:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1442:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/geometry"));
    1443:   PetscCall(VecView(newcoords, viewer));
    1444:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1445:   PetscCall(VecDestroy(&newcoords));
    1446:   PetscFunctionReturn(PETSC_SUCCESS);
    1447: }

    1449: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
    1450: {
    1451:   const char          *topologydm_name;
    1452:   const PetscInt      *gpoint;
    1453:   PetscInt             numLabels;
    1454:   PetscBool            omitCelltypes = PETSC_FALSE, ocompress;
    1455:   DMPlexStorageVersion version;
    1456:   char                 group[PETSC_MAX_PATH_LEN];

    1458:   PetscFunctionBegin;
    1459:   PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_omit_celltypes", &omitCelltypes, NULL));
    1460:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    1461:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
    1462:     PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
    1463:     PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
    1464:   }
    1465:   PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
    1466:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1467:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
    1468:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
    1469:   } else {
    1470:     PetscCall(PetscStrncpy(group, "/labels", sizeof(group)));
    1471:   }
    1472:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));
    1473:   PetscCall(DMGetNumLabels(dm, &numLabels));
    1474:   for (PetscInt l = 0; l < numLabels; ++l) {
    1475:     DMLabel         label;
    1476:     const char     *name;
    1477:     IS              valueIS, pvalueIS, globalValueIS;
    1478:     const PetscInt *values;
    1479:     PetscInt        numValues, v;
    1480:     PetscBool       isDepth, isCelltype, output;

    1482:     PetscCall(DMGetLabelByNum(dm, l, &label));
    1483:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
    1484:     PetscCall(DMGetLabelOutput(dm, name, &output));
    1485:     PetscCall(PetscStrncmp(name, "depth", 10, &isDepth));
    1486:     PetscCall(PetscStrncmp(name, "celltype", 10, &isCelltype));
    1487:     // TODO Should only filter out celltype if it can be calculated
    1488:     if (isDepth || (isCelltype && omitCelltypes) || !output) continue;
    1489:     PetscCall(PetscViewerHDF5PushGroup(viewer, name));
    1490:     PetscCall(DMLabelGetValueIS(label, &valueIS));
    1491:     /* Must copy to a new IS on the global comm */
    1492:     PetscCall(ISGetLocalSize(valueIS, &numValues));
    1493:     PetscCall(ISGetIndices(valueIS, &values));
    1494:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS));
    1495:     PetscCall(ISRestoreIndices(valueIS, &values));
    1496:     PetscCall(ISAllGather(pvalueIS, &globalValueIS));
    1497:     PetscCall(ISDestroy(&pvalueIS));
    1498:     PetscCall(ISSortRemoveDups(globalValueIS));
    1499:     PetscCall(ISGetLocalSize(globalValueIS, &numValues));
    1500:     PetscCall(ISGetIndices(globalValueIS, &values));
    1501:     for (v = 0; v < numValues; ++v) {
    1502:       IS              stratumIS, globalStratumIS;
    1503:       const PetscInt *spoints = NULL;
    1504:       PetscInt       *gspoints, n = 0, gn, p;
    1505:       const char     *iname = "indices";
    1506:       char            group[PETSC_MAX_PATH_LEN];

    1508:       PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v]));
    1509:       PetscCall(PetscViewerHDF5PushGroup(viewer, group));
    1510:       PetscCall(DMLabelGetStratumIS(label, values[v], &stratumIS));

    1512:       if (stratumIS) PetscCall(ISGetLocalSize(stratumIS, &n));
    1513:       if (stratumIS) PetscCall(ISGetIndices(stratumIS, &spoints));
    1514:       for (gn = 0, p = 0; p < n; ++p)
    1515:         if (gpoint[spoints[p]] >= 0) ++gn;
    1516:       PetscCall(PetscMalloc1(gn, &gspoints));
    1517:       for (gn = 0, p = 0; p < n; ++p)
    1518:         if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
    1519:       if (stratumIS) PetscCall(ISRestoreIndices(stratumIS, &spoints));
    1520:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS));
    1521:       PetscCall(PetscObjectSetName((PetscObject)globalStratumIS, iname));
    1522:       PetscCall(ISView(globalStratumIS, viewer));
    1523:       PetscCall(ISDestroy(&globalStratumIS));
    1524:       PetscCall(ISDestroy(&stratumIS));
    1525:       PetscCall(PetscViewerHDF5PopGroup(viewer));
    1526:     }
    1527:     PetscCall(ISRestoreIndices(globalValueIS, &values));
    1528:     PetscCall(ISDestroy(&globalValueIS));
    1529:     PetscCall(ISDestroy(&valueIS));
    1530:     PetscCall(PetscViewerHDF5PopGroup(viewer));
    1531:   }
    1532:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
    1533:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1534:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
    1535:   PetscFunctionReturn(PETSC_SUCCESS);
    1536: }

    1538: /* We only write cells and vertices. Does this screw up parallel reading? */
    1539: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
    1540: {
    1541:   IS                globalPointNumbers;
    1542:   PetscViewerFormat format;
    1543:   PetscBool         viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE, view_rank = PETSC_FALSE;

    1545:   PetscFunctionBegin;
    1546:   PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
    1547:   PetscCall(DMPlexCoordinatesView_HDF5_Internal(dm, viewer));

    1549:   PetscCall(PetscViewerGetFormat(viewer, &format));
    1550:   switch (format) {
    1551:   case PETSC_VIEWER_HDF5_VIZ:
    1552:     viz_geom  = PETSC_TRUE;
    1553:     xdmf_topo = PETSC_TRUE;
    1554:     view_rank = PETSC_TRUE;
    1555:     break;
    1556:   case PETSC_VIEWER_HDF5_XDMF:
    1557:     xdmf_topo = PETSC_TRUE;
    1558:     break;
    1559:   case PETSC_VIEWER_HDF5_PETSC:
    1560:     petsc_topo = PETSC_TRUE;
    1561:     break;
    1562:   case PETSC_VIEWER_DEFAULT:
    1563:   case PETSC_VIEWER_NATIVE:
    1564:     viz_geom   = PETSC_TRUE;
    1565:     xdmf_topo  = PETSC_TRUE;
    1566:     petsc_topo = PETSC_TRUE;
    1567:     break;
    1568:   default:
    1569:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
    1570:   }

    1572:   if (viz_geom) PetscCall(DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer));
    1573:   if (xdmf_topo) PetscCall(DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer));
    1574:   if (petsc_topo) {
    1575:     PetscBool viewLabels = PETSC_TRUE;

    1577:     PetscCall(DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer));
    1578:     PetscCall(PetscOptionsGetBool(NULL, dm->hdr.prefix, "-dm_plex_view_labels", &viewLabels, NULL));
    1579:     if (viewLabels) PetscCall(DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer));
    1580:   }
    1581:   if (view_rank) {
    1582:     Vec v;

    1584:     PetscCall(DMPlexCreateRankField(dm, &v));
    1585:     PetscCall(VecView(v, viewer));
    1586:     PetscCall(VecDestroy(&v));
    1587:   }
    1588:   PetscCall(ISDestroy(&globalPointNumbers));
    1589:   PetscFunctionReturn(PETSC_SUCCESS);
    1590: }

    1592: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
    1593: {
    1594:   DMPlexStorageVersion version;
    1595:   MPI_Comm             comm;
    1596:   const char          *topologydm_name;
    1597:   const char          *sectiondm_name;
    1598:   PetscSection         gsection;
    1599:   PetscBool            ocompress;

    1601:   PetscFunctionBegin;
    1602:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
    1603:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) {
    1604:     PetscCall(PetscViewerHDF5GetCompress(viewer, &ocompress));
    1605:     PetscCall(PetscViewerHDF5SetCompress(viewer, PETSC_TRUE));
    1606:   }
    1607:   PetscCall(PetscObjectGetComm((PetscObject)sectiondm, &comm));
    1608:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1609:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
    1610:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    1611:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    1612:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
    1613:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
    1614:   PetscCall(DMGetGlobalSection(sectiondm, &gsection));
    1615:   /* Save raw section */
    1616:   PetscCall(PetscSectionView(gsection, viewer));
    1617:   /* Save plex wrapper */
    1618:   {
    1619:     PetscInt        pStart, pEnd, p, n;
    1620:     IS              globalPointNumbers;
    1621:     const PetscInt *gpoints;
    1622:     IS              orderIS;
    1623:     PetscInt       *order;

    1625:     PetscCall(PetscSectionGetChart(gsection, &pStart, &pEnd));
    1626:     PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
    1627:     PetscCall(ISGetIndices(globalPointNumbers, &gpoints));
    1628:     for (p = pStart, n = 0; p < pEnd; ++p)
    1629:       if (gpoints[p] >= 0) n++;
    1630:     /* "order" is an array of global point numbers.
    1631:        When loading, it is used with topology/order array
    1632:        to match section points with plex topology points. */
    1633:     PetscCall(PetscMalloc1(n, &order));
    1634:     for (p = pStart, n = 0; p < pEnd; ++p)
    1635:       if (gpoints[p] >= 0) order[n++] = gpoints[p];
    1636:     PetscCall(ISRestoreIndices(globalPointNumbers, &gpoints));
    1637:     PetscCall(ISDestroy(&globalPointNumbers));
    1638:     PetscCall(ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS));
    1639:     PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
    1640:     PetscCall(ISView(orderIS, viewer));
    1641:     PetscCall(ISDestroy(&orderIS));
    1642:   }
    1643:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1644:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1645:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1646:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1647:   if (DMPlexStorageVersionGE(version, 3, 1, 0)) PetscCall(PetscViewerHDF5SetCompress(viewer, ocompress));
    1648:   PetscFunctionReturn(PETSC_SUCCESS);
    1649: }

    1651: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
    1652: {
    1653:   const char *topologydm_name;
    1654:   const char *sectiondm_name;
    1655:   const char *vec_name;
    1656:   PetscInt    bs;

    1658:   PetscFunctionBegin;
    1659:   /* Check consistency */
    1660:   {
    1661:     PetscSF pointsf, pointsf1;

    1663:     PetscCall(DMGetPointSF(dm, &pointsf));
    1664:     PetscCall(DMGetPointSF(sectiondm, &pointsf1));
    1665:     PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
    1666:   }
    1667:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1668:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
    1669:   PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
    1670:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    1671:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    1672:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
    1673:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
    1674:   PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
    1675:   PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
    1676:   PetscCall(VecGetBlockSize(vec, &bs));
    1677:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
    1678:   PetscCall(VecSetBlockSize(vec, 1));
    1679:   /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but,    */
    1680:   /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view    */
    1681:   /* is set to VecView_Plex, which would save vec in a predefined location.  */
    1682:   /* To save vec in where we want, we create a new Vec (temp) with           */
    1683:   /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
    1684:   {
    1685:     Vec                temp;
    1686:     const PetscScalar *array;
    1687:     PetscLayout        map;

    1689:     PetscCall(VecCreate(PetscObjectComm((PetscObject)vec), &temp));
    1690:     PetscCall(PetscObjectSetName((PetscObject)temp, vec_name));
    1691:     PetscCall(VecGetLayout(vec, &map));
    1692:     PetscCall(VecSetLayout(temp, map));
    1693:     PetscCall(VecSetUp(temp));
    1694:     PetscCall(VecGetArrayRead(vec, &array));
    1695:     PetscCall(VecPlaceArray(temp, array));
    1696:     PetscCall(VecView(temp, viewer));
    1697:     PetscCall(VecResetArray(temp));
    1698:     PetscCall(VecRestoreArrayRead(vec, &array));
    1699:     PetscCall(VecDestroy(&temp));
    1700:   }
    1701:   PetscCall(VecSetBlockSize(vec, bs));
    1702:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1703:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1704:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1705:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1706:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1707:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1708:   PetscFunctionReturn(PETSC_SUCCESS);
    1709: }

    1711: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
    1712: {
    1713:   MPI_Comm     comm;
    1714:   const char  *topologydm_name;
    1715:   const char  *sectiondm_name;
    1716:   const char  *vec_name;
    1717:   PetscSection section;
    1718:   PetscBool    includesConstraints;
    1719:   Vec          gvec;
    1720:   PetscInt     m, bs;

    1722:   PetscFunctionBegin;
    1723:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1724:   /* Check consistency */
    1725:   {
    1726:     PetscSF pointsf, pointsf1;

    1728:     PetscCall(DMGetPointSF(dm, &pointsf));
    1729:     PetscCall(DMGetPointSF(sectiondm, &pointsf1));
    1730:     PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
    1731:   }
    1732:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1733:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
    1734:   PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
    1735:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    1736:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    1737:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
    1738:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
    1739:   PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
    1740:   PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
    1741:   PetscCall(VecGetBlockSize(vec, &bs));
    1742:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
    1743:   PetscCall(VecCreate(comm, &gvec));
    1744:   PetscCall(PetscObjectSetName((PetscObject)gvec, vec_name));
    1745:   PetscCall(DMGetGlobalSection(sectiondm, &section));
    1746:   PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
    1747:   if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
    1748:   else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
    1749:   PetscCall(VecSetSizes(gvec, m, PETSC_DECIDE));
    1750:   PetscCall(VecSetUp(gvec));
    1751:   PetscCall(DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec));
    1752:   PetscCall(DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec));
    1753:   PetscCall(VecView(gvec, viewer));
    1754:   PetscCall(VecDestroy(&gvec));
    1755:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1756:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1757:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1758:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1759:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1760:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1761:   PetscFunctionReturn(PETSC_SUCCESS);
    1762: }

    1764: struct _n_LoadLabelsCtx {
    1765:   MPI_Comm    comm;
    1766:   PetscMPIInt rank;
    1767:   DM          dm;
    1768:   PetscViewer viewer;
    1769:   DMLabel     label;
    1770:   PetscSF     sfXC;
    1771:   PetscLayout layoutX;
    1772: };
    1773: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;

    1775: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
    1776: {
    1777:   PetscFunctionBegin;
    1778:   PetscCall(PetscNew(ctx));
    1779:   PetscCall(PetscObjectReference((PetscObject)((*ctx)->dm = dm)));
    1780:   PetscCall(PetscObjectReference((PetscObject)((*ctx)->viewer = viewer)));
    1781:   PetscCall(PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm));
    1782:   PetscCallMPI(MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank));
    1783:   (*ctx)->sfXC = sfXC;
    1784:   if (sfXC) {
    1785:     PetscInt nX;

    1787:     PetscCall(PetscObjectReference((PetscObject)sfXC));
    1788:     PetscCall(PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL));
    1789:     PetscCall(PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX));
    1790:   }
    1791:   PetscFunctionReturn(PETSC_SUCCESS);
    1792: }

    1794: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
    1795: {
    1796:   PetscFunctionBegin;
    1797:   if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
    1798:   PetscCall(DMDestroy(&(*ctx)->dm));
    1799:   PetscCall(PetscViewerDestroy(&(*ctx)->viewer));
    1800:   PetscCall(PetscSFDestroy(&(*ctx)->sfXC));
    1801:   PetscCall(PetscLayoutDestroy(&(*ctx)->layoutX));
    1802:   PetscCall(PetscFree(*ctx));
    1803:   PetscFunctionReturn(PETSC_SUCCESS);
    1804: }

    1806: /*
    1807:     A: on-disk points
    1808:     X: global points [0, NX)
    1809:     C: distributed plex points
    1810: */
    1811: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
    1812: {
    1813:   MPI_Comm        comm    = ctx->comm;
    1814:   PetscSF         sfXC    = ctx->sfXC;
    1815:   PetscLayout     layoutX = ctx->layoutX;
    1816:   PetscSF         sfXA;
    1817:   const PetscInt *A_points;
    1818:   PetscInt        nX, nC;
    1819:   PetscInt        n;

    1821:   PetscFunctionBegin;
    1822:   PetscCall(PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL));
    1823:   PetscCall(ISGetLocalSize(stratumIS, &n));
    1824:   PetscCall(ISGetIndices(stratumIS, &A_points));
    1825:   PetscCall(PetscSFCreate(comm, &sfXA));
    1826:   PetscCall(PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points));
    1827:   PetscCall(ISCreate(comm, newStratumIS));
    1828:   PetscCall(ISSetType(*newStratumIS, ISGENERAL));
    1829:   {
    1830:     PetscInt   i;
    1831:     PetscBool *A_mask, *X_mask, *C_mask;

    1833:     PetscCall(PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask));
    1834:     for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE;
    1835:     PetscCall(PetscSFReduceBegin(sfXA, MPI_C_BOOL, A_mask, X_mask, MPI_REPLACE));
    1836:     PetscCall(PetscSFReduceEnd(sfXA, MPI_C_BOOL, A_mask, X_mask, MPI_REPLACE));
    1837:     PetscCall(PetscSFBcastBegin(sfXC, MPI_C_BOOL, X_mask, C_mask, MPI_LOR));
    1838:     PetscCall(PetscSFBcastEnd(sfXC, MPI_C_BOOL, X_mask, C_mask, MPI_LOR));
    1839:     PetscCall(ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask));
    1840:     PetscCall(PetscFree3(A_mask, X_mask, C_mask));
    1841:   }
    1842:   PetscCall(PetscSFDestroy(&sfXA));
    1843:   PetscCall(ISRestoreIndices(stratumIS, &A_points));
    1844:   PetscFunctionReturn(PETSC_SUCCESS);
    1845: }

    1847: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
    1848: {
    1849:   LoadLabelsCtx   ctx    = (LoadLabelsCtx)op_data;
    1850:   PetscViewer     viewer = ctx->viewer;
    1851:   DMLabel         label  = ctx->label;
    1852:   MPI_Comm        comm   = ctx->comm;
    1853:   IS              stratumIS;
    1854:   const PetscInt *ind;
    1855:   PetscInt        value, N, i;

    1857:   PetscCall(PetscOptionsStringToInt(vname, &value));
    1858:   PetscCall(ISCreate(comm, &stratumIS));
    1859:   PetscCall(PetscObjectSetName((PetscObject)stratumIS, "indices"));
    1860:   PetscCall(PetscViewerHDF5PushGroup(viewer, vname)); /* labels/<lname>/<vname> */

    1862:   if (!ctx->sfXC) {
    1863:     /* Force serial load */
    1864:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N));
    1865:     PetscCall(PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0));
    1866:     PetscCall(PetscLayoutSetSize(stratumIS->map, N));
    1867:   }
    1868:   PetscCall(ISLoad(stratumIS, viewer));

    1870:   if (ctx->sfXC) {
    1871:     IS newStratumIS;

    1873:     PetscCallHDF5(ReadLabelStratumHDF5_Distribute_Private, (stratumIS, ctx, &newStratumIS));
    1874:     PetscCall(ISDestroy(&stratumIS));
    1875:     stratumIS = newStratumIS;
    1876:   }

    1878:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1879:   PetscCall(ISGetLocalSize(stratumIS, &N));
    1880:   PetscCall(ISGetIndices(stratumIS, &ind));
    1881:   for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, ind[i], value));
    1882:   PetscCall(ISRestoreIndices(stratumIS, &ind));
    1883:   PetscCall(ISDestroy(&stratumIS));
    1884:   return 0;
    1885: }

    1887: /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */
    1888: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
    1889: {
    1890:   LoadLabelsCtx  ctx = (LoadLabelsCtx)op_data;
    1891:   DM             dm  = ctx->dm;
    1892:   hsize_t        idx = 0;
    1893:   PetscErrorCode ierr;
    1894:   PetscBool      flg;
    1895:   herr_t         err;

    1897:   PetscCall(DMHasLabel(dm, lname, &flg));
    1898:   if (flg) PetscCall(DMRemoveLabel(dm, lname, NULL));
    1899:   ierr = DMCreateLabel(dm, lname);
    1900:   if (ierr) return (herr_t)ierr;
    1901:   ierr = DMGetLabel(dm, lname, &ctx->label);
    1902:   if (ierr) return (herr_t)ierr;
    1903:   ierr = PetscViewerHDF5PushGroup(ctx->viewer, lname);
    1904:   if (ierr) return (herr_t)ierr;
    1905:   /* Iterate over the label's strata */
    1906:   PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
    1907:   ierr = PetscViewerHDF5PopGroup(ctx->viewer);
    1908:   if (ierr) return (herr_t)ierr;
    1909:   return err;
    1910: }

    1912: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
    1913: {
    1914:   const char          *topologydm_name;
    1915:   LoadLabelsCtx        ctx;
    1916:   hsize_t              idx = 0;
    1917:   char                 group[PETSC_MAX_PATH_LEN];
    1918:   DMPlexStorageVersion version;
    1919:   PetscBool            distributed, hasGroup;

    1921:   PetscFunctionBegin;
    1922:   PetscCall(DMPlexIsDistributed(dm, &distributed));
    1923:   if (distributed) PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
    1924:   PetscCall(LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx));
    1925:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    1926:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
    1927:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
    1928:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
    1929:   } else {
    1930:     PetscCall(PetscStrncpy(group, "labels", sizeof(group)));
    1931:   }
    1932:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));
    1933:   PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup));
    1934:   if (hasGroup) {
    1935:     hid_t fileId, groupId;

    1937:     PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId));
    1938:     /* Iterate over labels */
    1939:     PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
    1940:     PetscCallHDF5(H5Gclose, (groupId));
    1941:   }
    1942:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    1943:   PetscCall(LoadLabelsCtxDestroy(&ctx));
    1944:   PetscFunctionReturn(PETSC_SUCCESS);
    1945: }

    1947: static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm)
    1948: {
    1949:   MPI_Comm        comm;
    1950:   PetscMPIInt     size, rank;
    1951:   PetscInt        dist_size;
    1952:   const char     *distribution_name;
    1953:   PetscInt        p, lsize;
    1954:   IS              chartSizesIS, ownersIS, gpointsIS;
    1955:   const PetscInt *chartSize, *owners, *gpoints;
    1956:   PetscLayout     layout;
    1957:   PetscBool       has;

    1959:   PetscFunctionBegin;
    1960:   *distsf = NULL;
    1961:   *distdm = NULL;
    1962:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1963:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1964:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1965:   PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
    1966:   if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
    1967:   PetscCall(PetscLogEventBegin(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
    1968:   PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &has));
    1969:   if (!has) {
    1970:     const char *full_group;

    1972:     PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &full_group));
    1973:     PetscCheck(has, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Distribution %s cannot be found: HDF5 group %s not found in file", distribution_name, full_group);
    1974:   }
    1975:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size));
    1976:   PetscCheck(dist_size == (PetscInt)size, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mismatching comm sizes: comm size of this session (%d) != comm size used for %s (%" PetscInt_FMT ")", size, distribution_name, dist_size);
    1977:   PetscCall(ISCreate(comm, &chartSizesIS));
    1978:   PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
    1979:   PetscCall(ISCreate(comm, &ownersIS));
    1980:   PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
    1981:   PetscCall(ISCreate(comm, &gpointsIS));
    1982:   PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
    1983:   PetscCall(PetscLayoutSetLocalSize(chartSizesIS->map, 1));
    1984:   PetscCall(ISLoad(chartSizesIS, viewer));
    1985:   PetscCall(ISGetIndices(chartSizesIS, &chartSize));
    1986:   PetscCall(PetscLayoutSetLocalSize(ownersIS->map, *chartSize));
    1987:   PetscCall(PetscLayoutSetLocalSize(gpointsIS->map, *chartSize));
    1988:   PetscCall(ISLoad(ownersIS, viewer));
    1989:   PetscCall(ISLoad(gpointsIS, viewer));
    1990:   PetscCall(ISGetIndices(ownersIS, &owners));
    1991:   PetscCall(ISGetIndices(gpointsIS, &gpoints));
    1992:   PetscCall(PetscSFCreate(comm, distsf));
    1993:   PetscCall(PetscSFSetFromOptions(*distsf));
    1994:   PetscCall(PetscLayoutCreate(comm, &layout));
    1995:   PetscCall(PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL));
    1996:   PetscCall(PetscLayoutSetLocalSize(layout, lsize));
    1997:   PetscCall(PetscLayoutSetBlockSize(layout, 1));
    1998:   PetscCall(PetscLayoutSetUp(layout));
    1999:   PetscCall(PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints));
    2000:   PetscCall(PetscLayoutDestroy(&layout));
    2001:   /* Migrate DM */
    2002:   {
    2003:     PetscInt     pStart, pEnd;
    2004:     PetscSFNode *buffer0, *buffer1, *buffer2;

    2006:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
    2007:     PetscCall(PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1));
    2008:     PetscCall(PetscMalloc1(*chartSize, &buffer2));
    2009:     {
    2010:       PetscSF            workPointSF;
    2011:       PetscInt           workNroots, workNleaves;
    2012:       const PetscInt    *workIlocal;
    2013:       const PetscSFNode *workIremote;

    2015:       for (p = pStart; p < pEnd; ++p) {
    2016:         buffer0[p - pStart].rank  = rank;
    2017:         buffer0[p - pStart].index = p - pStart;
    2018:       }
    2019:       PetscCall(DMGetPointSF(dm, &workPointSF));
    2020:       PetscCall(PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote));
    2021:       for (p = 0; p < workNleaves; ++p) {
    2022:         PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p);

    2024:         buffer0[workIlocalp].rank = -1;
    2025:       }
    2026:     }
    2027:     for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
    2028:     for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1;
    2029:     PetscCall(PetscSFReduceBegin(sf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
    2030:     PetscCall(PetscSFReduceEnd(sf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
    2031:     PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer2, MPI_REPLACE));
    2032:     PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer2, MPI_REPLACE));
    2033:     if (PetscDefined(USE_DEBUG)) {
    2034:       for (p = 0; p < *chartSize; ++p) PetscCheck(buffer2[p].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Found negative root rank %" PetscInt_FMT " at local point %" PetscInt_FMT " on rank %d when making migrationSF", buffer2[p].rank, p, rank);
    2035:     }
    2036:     PetscCall(PetscFree2(buffer0, buffer1));
    2037:     PetscCall(DMCreate(comm, distdm));
    2038:     PetscCall(DMSetType(*distdm, DMPLEX));
    2039:     PetscCall(PetscObjectSetName((PetscObject)*distdm, ((PetscObject)dm)->name));
    2040:     PetscCall(DMPlexDistributionSetName(*distdm, distribution_name));
    2041:     {
    2042:       PetscSF migrationSF;

    2044:       PetscCall(PetscSFCreate(comm, &migrationSF));
    2045:       PetscCall(PetscSFSetFromOptions(migrationSF));
    2046:       PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER));
    2047:       PetscCall(PetscSFSetUp(migrationSF));
    2048:       PetscCall(DMPlexMigrate(dm, migrationSF, *distdm));
    2049:       PetscCall(PetscSFDestroy(&migrationSF));
    2050:     }
    2051:   }
    2052:   /* Set pointSF */
    2053:   {
    2054:     PetscSF      pointSF;
    2055:     PetscInt    *ilocal, nleaves, q;
    2056:     PetscSFNode *iremote, *buffer0, *buffer1;

    2058:     PetscCall(PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1));
    2059:     for (p = 0, nleaves = 0; p < *chartSize; ++p) {
    2060:       if (owners[p] == rank) {
    2061:         buffer0[p].rank = rank;
    2062:       } else {
    2063:         buffer0[p].rank = -1;
    2064:         nleaves++;
    2065:       }
    2066:       buffer0[p].index = p;
    2067:     }
    2068:     for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
    2069:     PetscCall(PetscSFReduceBegin(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
    2070:     PetscCall(PetscSFReduceEnd(*distsf, MPIU_SF_NODE, buffer0, buffer1, MPI_MAXLOC));
    2071:     for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1;
    2072:     PetscCall(PetscSFBcastBegin(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE));
    2073:     PetscCall(PetscSFBcastEnd(*distsf, MPIU_SF_NODE, buffer1, buffer0, MPI_REPLACE));
    2074:     if (PetscDefined(USE_DEBUG)) {
    2075:       for (p = 0; p < *chartSize; ++p) PetscCheck(buffer0[p].rank >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Found negative root rank %" PetscInt_FMT " at local point %" PetscInt_FMT " on rank %d when making pointSF", buffer0[p].rank, p, rank);
    2076:     }
    2077:     PetscCall(PetscMalloc1(nleaves, &ilocal));
    2078:     PetscCall(PetscMalloc1(nleaves, &iremote));
    2079:     for (p = 0, q = 0; p < *chartSize; ++p) {
    2080:       if (buffer0[p].rank != rank) {
    2081:         ilocal[q]        = p;
    2082:         iremote[q].rank  = buffer0[p].rank;
    2083:         iremote[q].index = buffer0[p].index;
    2084:         q++;
    2085:       }
    2086:     }
    2087:     PetscCheck(q == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching leaf sizes: %" PetscInt_FMT " != %" PetscInt_FMT, q, nleaves);
    2088:     PetscCall(PetscFree2(buffer0, buffer1));
    2089:     PetscCall(PetscSFCreate(comm, &pointSF));
    2090:     PetscCall(PetscSFSetFromOptions(pointSF));
    2091:     PetscCall(PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
    2092:     PetscCall(DMSetPointSF(*distdm, pointSF));
    2093:     {
    2094:       DM cdm;

    2096:       PetscCall(DMGetCoordinateDM(*distdm, &cdm));
    2097:       PetscCall(DMSetPointSF(cdm, pointSF));
    2098:     }
    2099:     PetscCall(PetscSFDestroy(&pointSF));
    2100:   }
    2101:   PetscCall(ISRestoreIndices(chartSizesIS, &chartSize));
    2102:   PetscCall(ISRestoreIndices(ownersIS, &owners));
    2103:   PetscCall(ISRestoreIndices(gpointsIS, &gpoints));
    2104:   PetscCall(ISDestroy(&chartSizesIS));
    2105:   PetscCall(ISDestroy(&ownersIS));
    2106:   PetscCall(ISDestroy(&gpointsIS));
    2107:   /* Record that overlap has been manually created.               */
    2108:   /* This is to pass `DMPlexCheckPointSF()`, which checks that    */
    2109:   /* pointSF does not contain cells in the leaves if overlap = 0. */
    2110:   PetscCall(DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL));
    2111:   PetscCall(DMPlexDistributeSetDefault(*distdm, PETSC_FALSE));
    2112:   PetscCall(DMPlexReorderSetDefault(*distdm, DM_REORDER_DEFAULT_FALSE));
    2113:   PetscCall(PetscLogEventEnd(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
    2114:   PetscFunctionReturn(PETSC_SUCCESS);
    2115: }

    2117: // Serial load of topology
    2118: static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf)
    2119: {
    2120:   MPI_Comm        comm;
    2121:   const char     *pointsName, *coneSizesName, *conesName, *orientationsName;
    2122:   IS              pointsIS, coneSizesIS, conesIS, orientationsIS;
    2123:   const PetscInt *points, *coneSizes, *cones, *orientations;
    2124:   PetscInt       *cone, *ornt;
    2125:   PetscInt        dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
    2126:   PetscMPIInt     size, rank;

    2128:   PetscFunctionBegin;
    2129:   pointsName       = "order";
    2130:   coneSizesName    = "cones";
    2131:   conesName        = "cells";
    2132:   orientationsName = "orientation";
    2133:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    2134:   PetscCallMPI(MPI_Comm_size(comm, &size));
    2135:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    2136:   PetscCall(ISCreate(comm, &pointsIS));
    2137:   PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
    2138:   PetscCall(ISCreate(comm, &coneSizesIS));
    2139:   PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
    2140:   PetscCall(ISCreate(comm, &conesIS));
    2141:   PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
    2142:   PetscCall(ISCreate(comm, &orientationsIS));
    2143:   PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
    2144:   PetscCall(PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim));
    2145:   PetscCall(DMSetDimension(dm, dim));
    2146:   {
    2147:     /* Force serial load */
    2148:     PetscCall(PetscInfo(dm, "Loading DM %s in serial\n", dm->hdr.name));
    2149:     PetscCall(PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np));
    2150:     PetscCall(PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0));
    2151:     PetscCall(PetscLayoutSetSize(pointsIS->map, Np));
    2152:     pEnd = rank == 0 ? Np : 0;
    2153:     PetscCall(PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np));
    2154:     PetscCall(PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0));
    2155:     PetscCall(PetscLayoutSetSize(coneSizesIS->map, Np));
    2156:     PetscCall(PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N));
    2157:     PetscCall(PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0));
    2158:     PetscCall(PetscLayoutSetSize(conesIS->map, N));
    2159:     PetscCall(PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N));
    2160:     PetscCall(PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0));
    2161:     PetscCall(PetscLayoutSetSize(orientationsIS->map, N));
    2162:   }
    2163:   PetscCall(ISLoad(pointsIS, viewer));
    2164:   PetscCall(ISLoad(coneSizesIS, viewer));
    2165:   PetscCall(ISLoad(conesIS, viewer));
    2166:   PetscCall(ISLoad(orientationsIS, viewer));
    2167:   /* Create Plex */
    2168:   PetscCall(DMPlexSetChart(dm, 0, pEnd));
    2169:   PetscCall(ISGetIndices(pointsIS, &points));
    2170:   PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
    2171:   for (p = 0; p < pEnd; ++p) {
    2172:     PetscCall(DMPlexSetConeSize(dm, points[p], coneSizes[p]));
    2173:     maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
    2174:   }
    2175:   PetscCall(DMSetUp(dm));
    2176:   PetscCall(ISGetIndices(conesIS, &cones));
    2177:   PetscCall(ISGetIndices(orientationsIS, &orientations));
    2178:   PetscCall(PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt));
    2179:   for (p = 0, q = 0; p < pEnd; ++p) {
    2180:     for (c = 0; c < coneSizes[p]; ++c, ++q) {
    2181:       cone[c] = cones[q];
    2182:       ornt[c] = orientations[q];
    2183:     }
    2184:     PetscCall(DMPlexSetCone(dm, points[p], cone));
    2185:     PetscCall(DMPlexSetConeOrientation(dm, points[p], ornt));
    2186:   }
    2187:   PetscCall(PetscFree2(cone, ornt));
    2188:   /* Create global section migration SF */
    2189:   if (sf) {
    2190:     PetscLayout layout;
    2191:     PetscInt   *globalIndices;

    2193:     PetscCall(PetscMalloc1(pEnd, &globalIndices));
    2194:     /* plex point == globalPointNumber in this case */
    2195:     for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
    2196:     PetscCall(PetscLayoutCreate(comm, &layout));
    2197:     PetscCall(PetscLayoutSetSize(layout, Np));
    2198:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
    2199:     PetscCall(PetscLayoutSetUp(layout));
    2200:     PetscCall(PetscSFCreate(comm, sf));
    2201:     PetscCall(PetscSFSetFromOptions(*sf));
    2202:     PetscCall(PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices));
    2203:     PetscCall(PetscLayoutDestroy(&layout));
    2204:     PetscCall(PetscFree(globalIndices));
    2205:   }
    2206:   /* Clean-up */
    2207:   PetscCall(ISRestoreIndices(pointsIS, &points));
    2208:   PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
    2209:   PetscCall(ISRestoreIndices(conesIS, &cones));
    2210:   PetscCall(ISRestoreIndices(orientationsIS, &orientations));
    2211:   PetscCall(ISDestroy(&pointsIS));
    2212:   PetscCall(ISDestroy(&coneSizesIS));
    2213:   PetscCall(ISDestroy(&conesIS));
    2214:   PetscCall(ISDestroy(&orientationsIS));
    2215:   /* Fill in the rest of the topology structure */
    2216:   PetscCall(DMPlexSymmetrize(dm));
    2217:   PetscCall(DMPlexStratify(dm));
    2218:   PetscFunctionReturn(PETSC_SUCCESS);
    2219: }

    2221: /* Representation of two DMPlex strata in 0-based global numbering */
    2222: struct _n_PlexLayer {
    2223:   PetscInt     d;
    2224:   IS           conesIS, orientationsIS;
    2225:   PetscSection coneSizesSection;
    2226:   PetscLayout  vertexLayout;
    2227:   PetscSF      overlapSF, l2gSF; //TODO maybe confusing names (in DMPlex in general)
    2228:   PetscInt     offset, conesOffset, leafOffset;
    2229: };
    2230: typedef struct _n_PlexLayer *PlexLayer;

    2232: static PetscErrorCode PlexLayerDestroy(PlexLayer *layer)
    2233: {
    2234:   PetscFunctionBegin;
    2235:   if (!*layer) PetscFunctionReturn(PETSC_SUCCESS);
    2236:   PetscCall(PetscSectionDestroy(&(*layer)->coneSizesSection));
    2237:   PetscCall(ISDestroy(&(*layer)->conesIS));
    2238:   PetscCall(ISDestroy(&(*layer)->orientationsIS));
    2239:   PetscCall(PetscSFDestroy(&(*layer)->overlapSF));
    2240:   PetscCall(PetscSFDestroy(&(*layer)->l2gSF));
    2241:   PetscCall(PetscLayoutDestroy(&(*layer)->vertexLayout));
    2242:   PetscCall(PetscFree(*layer));
    2243:   PetscFunctionReturn(PETSC_SUCCESS);
    2244: }

    2246: static PetscErrorCode PlexLayerCreate_Private(PlexLayer *layer)
    2247: {
    2248:   PetscFunctionBegin;
    2249:   PetscCall(PetscNew(layer));
    2250:   (*layer)->d           = -1;
    2251:   (*layer)->offset      = -1;
    2252:   (*layer)->conesOffset = -1;
    2253:   (*layer)->leafOffset  = -1;
    2254:   PetscFunctionReturn(PETSC_SUCCESS);
    2255: }

    2257: // Parallel load of a depth stratum
    2258: static PetscErrorCode PlexLayerLoad_Private(PlexLayer layer, PetscViewer viewer, PetscInt d, PetscLayout pointsLayout)
    2259: {
    2260:   char         path[128];
    2261:   MPI_Comm     comm;
    2262:   const char  *coneSizesName, *conesName, *orientationsName;
    2263:   IS           coneSizesIS, conesIS, orientationsIS;
    2264:   PetscSection coneSizesSection;
    2265:   PetscLayout  vertexLayout = NULL;
    2266:   PetscInt     s;

    2268:   PetscFunctionBegin;
    2269:   coneSizesName    = "cone_sizes";
    2270:   conesName        = "cones";
    2271:   orientationsName = "orientations";
    2272:   PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));

    2274:   /* query size of next lower depth stratum (next lower dimension) */
    2275:   if (d > 0) {
    2276:     PetscInt NVertices;

    2278:     PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT "/%s", d - 1, coneSizesName));
    2279:     PetscCall(PetscViewerHDF5ReadSizes(viewer, path, NULL, &NVertices));
    2280:     PetscCall(PetscLayoutCreate(comm, &vertexLayout));
    2281:     PetscCall(PetscLayoutSetSize(vertexLayout, NVertices));
    2282:     PetscCall(PetscLayoutSetUp(vertexLayout));
    2283:   }

    2285:   PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT, d));
    2286:   PetscCall(PetscViewerHDF5PushGroup(viewer, path));

    2288:   /* create coneSizesSection from stored IS coneSizes */
    2289:   {
    2290:     const PetscInt *coneSizes;

    2292:     PetscCall(ISCreate(comm, &coneSizesIS));
    2293:     PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
    2294:     if (pointsLayout) PetscCall(ISSetLayout(coneSizesIS, pointsLayout));
    2295:     PetscCall(ISLoad(coneSizesIS, viewer));
    2296:     if (!pointsLayout) PetscCall(ISGetLayout(coneSizesIS, &pointsLayout));
    2297:     PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
    2298:     PetscCall(PetscSectionCreate(comm, &coneSizesSection));
    2299:     //TODO different start ?
    2300:     PetscCall(PetscSectionSetChart(coneSizesSection, 0, pointsLayout->n));
    2301:     for (s = 0; s < pointsLayout->n; ++s) PetscCall(PetscSectionSetDof(coneSizesSection, s, coneSizes[s]));
    2302:     PetscCall(PetscSectionSetUp(coneSizesSection));
    2303:     PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
    2304:     {
    2305:       PetscLayout tmp = NULL;

    2307:       /* We need to keep the layout until the end of function */
    2308:       PetscCall(PetscLayoutReference((PetscLayout)pointsLayout, &tmp));
    2309:     }
    2310:     PetscCall(ISDestroy(&coneSizesIS));
    2311:   }

    2313:   /* use value layout of coneSizesSection as layout of cones and orientations */
    2314:   {
    2315:     PetscLayout conesLayout;

    2317:     PetscCall(PetscSectionGetValueLayout(comm, coneSizesSection, &conesLayout));
    2318:     PetscCall(ISCreate(comm, &conesIS));
    2319:     PetscCall(ISCreate(comm, &orientationsIS));
    2320:     PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
    2321:     PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
    2322:     PetscCall(PetscLayoutDuplicate(conesLayout, &conesIS->map));
    2323:     PetscCall(PetscLayoutDuplicate(conesLayout, &orientationsIS->map));
    2324:     PetscCall(ISLoad(conesIS, viewer));
    2325:     PetscCall(ISLoad(orientationsIS, viewer));
    2326:     PetscCall(PetscLayoutDestroy(&conesLayout));
    2327:   }

    2329:   /* check assertion that layout of points is the same as point layout of coneSizesSection */
    2330:   {
    2331:     PetscLayout pointsLayout0;
    2332:     PetscBool   flg;

    2334:     PetscCall(PetscSectionGetPointLayout(comm, coneSizesSection, &pointsLayout0));
    2335:     PetscCall(PetscLayoutCompare(pointsLayout, pointsLayout0, &flg));
    2336:     PetscCheck(flg, comm, PETSC_ERR_PLIB, "points layout != coneSizesSection point layout");
    2337:     PetscCall(PetscLayoutDestroy(&pointsLayout0));
    2338:   }
    2339:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    2340:   PetscCall(PetscLayoutDestroy(&pointsLayout));

    2342:   layer->d                = d;
    2343:   layer->conesIS          = conesIS;
    2344:   layer->coneSizesSection = coneSizesSection;
    2345:   layer->orientationsIS   = orientationsIS;
    2346:   layer->vertexLayout     = vertexLayout;
    2347:   PetscFunctionReturn(PETSC_SUCCESS);
    2348: }

    2350: static PetscErrorCode PlexLayerDistribute_Private(PlexLayer layer, PetscSF cellLocalToGlobalSF)
    2351: {
    2352:   IS           newConesIS, newOrientationsIS;
    2353:   PetscSection newConeSizesSection;
    2354:   MPI_Comm     comm;

    2356:   PetscFunctionBegin;
    2357:   PetscCall(PetscObjectGetComm((PetscObject)cellLocalToGlobalSF, &comm));
    2358:   PetscCall(PetscSectionCreate(comm, &newConeSizesSection));
    2359:   //TODO rename to something like ISDistribute() with optional PetscSection argument
    2360:   PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->conesIS, newConeSizesSection, &newConesIS));
    2361:   PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->orientationsIS, newConeSizesSection, &newOrientationsIS));

    2363:   PetscCall(PetscObjectSetName((PetscObject)newConeSizesSection, ((PetscObject)layer->coneSizesSection)->name));
    2364:   PetscCall(PetscObjectSetName((PetscObject)newConesIS, ((PetscObject)layer->conesIS)->name));
    2365:   PetscCall(PetscObjectSetName((PetscObject)newOrientationsIS, ((PetscObject)layer->orientationsIS)->name));
    2366:   PetscCall(PetscSectionDestroy(&layer->coneSizesSection));
    2367:   PetscCall(ISDestroy(&layer->conesIS));
    2368:   PetscCall(ISDestroy(&layer->orientationsIS));
    2369:   layer->coneSizesSection = newConeSizesSection;
    2370:   layer->conesIS          = newConesIS;
    2371:   layer->orientationsIS   = newOrientationsIS;
    2372:   PetscFunctionReturn(PETSC_SUCCESS);
    2373: }

    2375: //TODO share code with DMPlexBuildFromCellListParallel()
    2376: #include <petsc/private/hashseti.h>
    2377: static PetscErrorCode PlexLayerCreateSFs_Private(PlexLayer layer, PetscSF *vertexOverlapSF, PetscSF *sfXC)
    2378: {
    2379:   PetscLayout     vertexLayout     = layer->vertexLayout;
    2380:   PetscSection    coneSection      = layer->coneSizesSection;
    2381:   IS              cellVertexData   = layer->conesIS;
    2382:   IS              coneOrientations = layer->orientationsIS;
    2383:   PetscSF         vl2gSF, vOverlapSF;
    2384:   PetscInt       *verticesAdj;
    2385:   PetscInt        i, n, numVerticesAdj;
    2386:   const PetscInt *cvd, *co = NULL;
    2387:   MPI_Comm        comm;

    2389:   PetscFunctionBegin;
    2390:   PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
    2391:   PetscCall(PetscSectionGetStorageSize(coneSection, &n));
    2392:   {
    2393:     PetscInt n0;

    2395:     PetscCall(ISGetLocalSize(cellVertexData, &n0));
    2396:     PetscCheck(n == n0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local size of IS cellVertexData = %" PetscInt_FMT " != %" PetscInt_FMT " = storage size of PetscSection coneSection", n0, n);
    2397:     PetscCall(ISGetIndices(cellVertexData, &cvd));
    2398:   }
    2399:   if (coneOrientations) {
    2400:     PetscInt n0;

    2402:     PetscCall(ISGetLocalSize(coneOrientations, &n0));
    2403:     PetscCheck(n == n0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local size of IS coneOrientations = %" PetscInt_FMT " != %" PetscInt_FMT " = storage size of PetscSection coneSection", n0, n);
    2404:     PetscCall(ISGetIndices(coneOrientations, &co));
    2405:   }
    2406:   /* Get/check global number of vertices */
    2407:   {
    2408:     PetscInt NVerticesInCells = PETSC_INT_MIN;

    2410:     /* NVerticesInCells = max(cellVertexData) + 1 */
    2411:     for (i = 0; i < n; i++)
    2412:       if (cvd[i] > NVerticesInCells) NVerticesInCells = cvd[i];
    2413:     ++NVerticesInCells;
    2414:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, comm));

    2416:     if (vertexLayout->n == PETSC_DECIDE && vertexLayout->N == PETSC_DECIDE) vertexLayout->N = NVerticesInCells;
    2417:     else
    2418:       PetscCheck(vertexLayout->N == PETSC_DECIDE || vertexLayout->N >= NVerticesInCells, comm, PETSC_ERR_ARG_SIZ, "Specified global number of vertices %" PetscInt_FMT " must be greater than or equal to the number of unique vertices in the cell-vertex dataset %" PetscInt_FMT,
    2419:                  vertexLayout->N, NVerticesInCells);
    2420:     PetscCall(PetscLayoutSetUp(vertexLayout));
    2421:   }
    2422:   /* Find locally unique vertices in cellVertexData */
    2423:   {
    2424:     PetscHSetI vhash;
    2425:     PetscInt   off = 0;

    2427:     PetscCall(PetscHSetICreate(&vhash));
    2428:     for (i = 0; i < n; ++i) PetscCall(PetscHSetIAdd(vhash, cvd[i]));
    2429:     PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj));
    2430:     PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj));
    2431:     PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj));
    2432:     PetscAssert(off == numVerticesAdj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assertion failed: off == numVerticesAdj (%" PetscInt_FMT " != %" PetscInt_FMT ")", off, numVerticesAdj);
    2433:     PetscCall(PetscHSetIDestroy(&vhash));
    2434:   }
    2435:   /* We must sort vertices to preserve numbering */
    2436:   PetscCall(PetscSortInt(numVerticesAdj, verticesAdj));
    2437:   /* Connect locally unique vertices with star forest */
    2438:   PetscCall(PetscSFCreateByMatchingIndices(vertexLayout, numVerticesAdj, verticesAdj, NULL, 0, numVerticesAdj, verticesAdj, NULL, 0, &vl2gSF, &vOverlapSF));
    2439:   PetscCall(PetscObjectSetName((PetscObject)vOverlapSF, "overlapSF"));
    2440:   PetscCall(PetscObjectSetName((PetscObject)vl2gSF, "localToGlobalSF"));

    2442:   PetscCall(PetscFree(verticesAdj));
    2443:   *vertexOverlapSF = vOverlapSF;
    2444:   *sfXC            = vl2gSF;
    2445:   PetscFunctionReturn(PETSC_SUCCESS);
    2446: }

    2448: static PetscErrorCode PlexLayerCreateCellSFs_Private(PlexLayer layer, PetscSF *cellOverlapSF, PetscSF *cellLocalToGlobalSF)
    2449: {
    2450:   PetscSection coneSection = layer->coneSizesSection;
    2451:   PetscInt     nCells;
    2452:   MPI_Comm     comm;

    2454:   PetscFunctionBegin;
    2455:   PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
    2456:   {
    2457:     PetscInt cStart;

    2459:     PetscCall(PetscSectionGetChart(coneSection, &cStart, &nCells));
    2460:     PetscCheck(cStart == 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "coneSection must start at 0");
    2461:   }
    2462:   /* Create overlapSF as empty SF with the right number of roots */
    2463:   PetscCall(PetscSFCreate(comm, cellOverlapSF));
    2464:   PetscCall(PetscSFSetGraph(*cellOverlapSF, nCells, 0, NULL, PETSC_USE_POINTER, NULL, PETSC_USE_POINTER));
    2465:   PetscCall(PetscSFSetUp(*cellOverlapSF));
    2466:   /* Create localToGlobalSF as identity mapping */
    2467:   {
    2468:     PetscLayout map;

    2470:     PetscCall(PetscLayoutCreateFromSizes(comm, nCells, PETSC_DECIDE, 1, &map));
    2471:     PetscCall(PetscSFCreateFromLayouts(map, map, cellLocalToGlobalSF));
    2472:     PetscCall(PetscSFSetUp(*cellLocalToGlobalSF));
    2473:     PetscCall(PetscLayoutDestroy(&map));
    2474:   }
    2475:   PetscFunctionReturn(PETSC_SUCCESS);
    2476: }

    2478: static PetscErrorCode DMPlexTopologyBuildFromLayers_Private(DM dm, PetscInt depth, PlexLayer *layers, IS strataPermutation)
    2479: {
    2480:   const PetscInt *permArr;
    2481:   PetscInt        d, nPoints;
    2482:   MPI_Comm        comm;

    2484:   PetscFunctionBegin;
    2485:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    2486:   PetscCall(ISGetIndices(strataPermutation, &permArr));

    2488:   /* Count points, strata offsets and cones offsets (taking strataPermutation into account) */
    2489:   {
    2490:     PetscInt stratumOffset = 0;
    2491:     PetscInt conesOffset   = 0;

    2493:     for (d = 0; d <= depth; d++) {
    2494:       const PetscInt  e = permArr[d];
    2495:       const PlexLayer l = layers[e];
    2496:       PetscInt        lo, n, size;

    2498:       PetscCall(PetscSectionGetChart(l->coneSizesSection, &lo, &n));
    2499:       PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &size));
    2500:       PetscCheck(lo == 0, comm, PETSC_ERR_PLIB, "starting point should be 0 in coneSizesSection %" PetscInt_FMT, d);
    2501:       l->offset      = stratumOffset;
    2502:       l->conesOffset = conesOffset;
    2503:       stratumOffset += n;
    2504:       conesOffset += size;
    2505:     }
    2506:     nPoints = stratumOffset;
    2507:   }

    2509:   /* Set interval for all plex points */
    2510:   //TODO we should store starting point of plex
    2511:   PetscCall(DMPlexSetChart(dm, 0, nPoints));

    2513:   /* Set up plex coneSection from layer coneSections */
    2514:   {
    2515:     PetscSection coneSection;

    2517:     PetscCall(DMPlexGetConeSection(dm, &coneSection));
    2518:     for (d = 0; d <= depth; d++) {
    2519:       const PlexLayer l = layers[d];
    2520:       PetscInt        n, q;

    2522:       PetscCall(PetscSectionGetChart(l->coneSizesSection, NULL, &n));
    2523:       for (q = 0; q < n; q++) {
    2524:         const PetscInt p = l->offset + q;
    2525:         PetscInt       coneSize;

    2527:         PetscCall(PetscSectionGetDof(l->coneSizesSection, q, &coneSize));
    2528:         PetscCall(PetscSectionSetDof(coneSection, p, coneSize));
    2529:       }
    2530:     }
    2531:   }
    2532:   //TODO this is terrible, DMSetUp_Plex() should be DMPlexSetUpSections() or so
    2533:   PetscCall(DMSetUp(dm));

    2535:   /* Renumber cones points from layer-global numbering to plex-local numbering */
    2536:   {
    2537:     PetscInt *cones, *ornts;

    2539:     PetscCall(DMPlexGetCones(dm, &cones));
    2540:     PetscCall(DMPlexGetConeOrientations(dm, &ornts));
    2541:     for (d = 1; d <= depth; d++) {
    2542:       const PlexLayer l = layers[d];
    2543:       PetscInt        i, lConesSize;
    2544:       PetscInt       *lCones;
    2545:       const PetscInt *lOrnts;
    2546:       PetscInt       *pCones = &cones[l->conesOffset];
    2547:       PetscInt       *pOrnts = &ornts[l->conesOffset];

    2549:       PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &lConesSize));
    2550:       /* Get cones in local plex numbering */
    2551:       {
    2552:         ISLocalToGlobalMapping l2g;
    2553:         PetscLayout            vertexLayout = l->vertexLayout;
    2554:         PetscSF                vertexSF     = layers[d - 1]->l2gSF; /* vertices of this layer are cells of previous layer */
    2555:         const PetscInt        *gCones;
    2556:         PetscInt               lConesSize0;

    2558:         PetscCall(ISGetLocalSize(l->conesIS, &lConesSize0));
    2559:         PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(conesIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);
    2560:         PetscCall(ISGetLocalSize(l->orientationsIS, &lConesSize0));
    2561:         PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(orientationsIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);

    2563:         PetscCall(PetscMalloc1(lConesSize, &lCones));
    2564:         PetscCall(ISGetIndices(l->conesIS, &gCones));
    2565:         PetscCall(ISLocalToGlobalMappingCreateSF(vertexSF, vertexLayout->rstart, &l2g));
    2566:         PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, lConesSize, gCones, &lConesSize0, lCones));
    2567:         PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "global to local does not cover all indices (%" PetscInt_FMT " of %" PetscInt_FMT ")", lConesSize0, lConesSize);
    2568:         PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
    2569:         PetscCall(ISRestoreIndices(l->conesIS, &gCones));
    2570:       }
    2571:       PetscCall(ISGetIndices(l->orientationsIS, &lOrnts));
    2572:       /* Set cones, need to add stratum offset */
    2573:       for (i = 0; i < lConesSize; i++) {
    2574:         pCones[i] = lCones[i] + layers[d - 1]->offset; /* cone points of current layer are points of previous layer */
    2575:         pOrnts[i] = lOrnts[i];
    2576:       }
    2577:       PetscCall(PetscFree(lCones));
    2578:       PetscCall(ISRestoreIndices(l->orientationsIS, &lOrnts));
    2579:     }
    2580:   }
    2581:   PetscCall(DMPlexSymmetrize(dm));
    2582:   PetscCall(DMPlexStratify(dm));
    2583:   PetscCall(ISRestoreIndices(strataPermutation, &permArr));
    2584:   PetscFunctionReturn(PETSC_SUCCESS);
    2585: }

    2587: static PetscErrorCode PlexLayerConcatenateSFs_Private(MPI_Comm comm, PetscInt depth, PlexLayer layers[], IS strataPermutation, PetscSF *overlapSF, PetscSF *l2gSF)
    2588: {
    2589:   PetscInt        d;
    2590:   PetscSF        *osfs, *lsfs;
    2591:   PetscInt       *leafOffsets;
    2592:   const PetscInt *permArr;

    2594:   PetscFunctionBegin;
    2595:   PetscCall(ISGetIndices(strataPermutation, &permArr));
    2596:   PetscCall(PetscCalloc3(depth + 1, &osfs, depth + 1, &lsfs, depth + 1, &leafOffsets));
    2597:   for (d = 0; d <= depth; d++) {
    2598:     const PetscInt e = permArr[d];

    2600:     PetscAssert(e == layers[e]->d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "assertion: e == layers[e]->d");
    2601:     osfs[d]        = layers[e]->overlapSF;
    2602:     lsfs[d]        = layers[e]->l2gSF;
    2603:     leafOffsets[d] = layers[e]->offset;
    2604:   }
    2605:   PetscCall(PetscSFConcatenate(comm, depth + 1, osfs, PETSCSF_CONCATENATE_ROOTMODE_LOCAL, leafOffsets, overlapSF));
    2606:   PetscCall(PetscSFConcatenate(comm, depth + 1, lsfs, PETSCSF_CONCATENATE_ROOTMODE_GLOBAL, leafOffsets, l2gSF));
    2607:   PetscCall(PetscFree3(osfs, lsfs, leafOffsets));
    2608:   PetscCall(ISRestoreIndices(strataPermutation, &permArr));
    2609:   PetscFunctionReturn(PETSC_SUCCESS);
    2610: }

    2612: // Parallel load of topology
    2613: static PetscErrorCode DMPlexTopologyLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF *sfXC)
    2614: {
    2615:   PlexLayer  *layers;
    2616:   IS          strataPermutation;
    2617:   PetscLayout pointsLayout;
    2618:   PetscInt    depth;
    2619:   PetscInt    d;
    2620:   MPI_Comm    comm;

    2622:   PetscFunctionBegin;
    2623:   {
    2624:     PetscInt dim;

    2626:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "depth", PETSC_INT, NULL, &depth));
    2627:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "cell_dim", PETSC_INT, NULL, &dim));
    2628:     PetscCall(DMSetDimension(dm, dim));
    2629:   }
    2630:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));

    2632:   PetscCall(PetscInfo(dm, "Loading DM %s in parallel\n", dm->hdr.name));
    2633:   {
    2634:     IS spOnComm;

    2636:     PetscCall(ISCreate(comm, &spOnComm));
    2637:     PetscCall(PetscObjectSetName((PetscObject)spOnComm, "permutation"));
    2638:     PetscCall(ISLoad(spOnComm, viewer));
    2639:     /* have the same serial IS on every rank */
    2640:     PetscCall(ISAllGather(spOnComm, &strataPermutation));
    2641:     PetscCall(PetscObjectSetName((PetscObject)strataPermutation, ((PetscObject)spOnComm)->name));
    2642:     PetscCall(ISDestroy(&spOnComm));
    2643:   }

    2645:   /* Create layers, load raw data for each layer */
    2646:   PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
    2647:   PetscCall(PetscMalloc1(depth + 1, &layers));
    2648:   for (d = depth, pointsLayout = NULL; d >= 0; pointsLayout = layers[d]->vertexLayout, d--) {
    2649:     PetscCall(PlexLayerCreate_Private(&layers[d]));
    2650:     PetscCall(PlexLayerLoad_Private(layers[d], viewer, d, pointsLayout));
    2651:   }
    2652:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */

    2654:   for (d = depth; d >= 0; d--) {
    2655:     /* Redistribute cells and vertices for each applicable layer */
    2656:     if (d < depth) PetscCall(PlexLayerDistribute_Private(layers[d], layers[d]->l2gSF));
    2657:     /* Create vertex overlap SF and vertex localToGlobal SF for each applicable layer */
    2658:     if (d > 0) PetscCall(PlexLayerCreateSFs_Private(layers[d], &layers[d - 1]->overlapSF, &layers[d - 1]->l2gSF));
    2659:   }
    2660:   /* Build trivial SFs for the cell layer as well */
    2661:   PetscCall(PlexLayerCreateCellSFs_Private(layers[depth], &layers[depth]->overlapSF, &layers[depth]->l2gSF));

    2663:   /* Build DMPlex topology from the layers */
    2664:   PetscCall(DMPlexTopologyBuildFromLayers_Private(dm, depth, layers, strataPermutation));

    2666:   /* Build overall point SF alias overlap SF */
    2667:   {
    2668:     PetscSF overlapSF;

    2670:     PetscCall(PlexLayerConcatenateSFs_Private(comm, depth, layers, strataPermutation, &overlapSF, sfXC));
    2671:     PetscCall(DMSetPointSF(dm, overlapSF));
    2672:     PetscCall(PetscSFDestroy(&overlapSF));
    2673:   }

    2675:   for (d = depth; d >= 0; d--) PetscCall(PlexLayerDestroy(&layers[d]));
    2676:   PetscCall(PetscFree(layers));
    2677:   PetscCall(ISDestroy(&strataPermutation));
    2678:   PetscFunctionReturn(PETSC_SUCCESS);
    2679: }

    2681: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC)
    2682: {
    2683:   DMPlexStorageVersion version;
    2684:   const char          *topologydm_name;
    2685:   char                 group[PETSC_MAX_PATH_LEN];
    2686:   PetscSF              sfwork = NULL;

    2688:   PetscFunctionBegin;
    2689:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    2690:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
    2691:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
    2692:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
    2693:   } else {
    2694:     PetscCall(PetscStrncpy(group, "/", sizeof(group)));
    2695:   }
    2696:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));

    2698:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
    2699:   if (version->major < 3) {
    2700:     PetscCall(DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork));
    2701:   } else {
    2702:     /* since DMPlexStorageVersion 3.0.0 */
    2703:     PetscCall(DMPlexTopologyLoad_HDF5_Private(dm, viewer, &sfwork));
    2704:   }
    2705:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */

    2707:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
    2708:     DM          distdm;
    2709:     PetscSF     distsf;
    2710:     const char *distribution_name;
    2711:     PetscBool   exists;

    2713:     PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
    2714:     /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */
    2715:     PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
    2716:     PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &exists));
    2717:     if (exists) {
    2718:       PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
    2719:       PetscCall(DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm));
    2720:       if (distdm) {
    2721:         PetscCall(DMPlexReplace_Internal(dm, &distdm));
    2722:         PetscCall(PetscSFDestroy(&sfwork));
    2723:         sfwork = distsf;
    2724:       }
    2725:       PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
    2726:     }
    2727:     PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
    2728:   }
    2729:   if (sfXC) {
    2730:     *sfXC = sfwork;
    2731:   } else {
    2732:     PetscCall(PetscSFDestroy(&sfwork));
    2733:   }

    2735:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    2736:   PetscFunctionReturn(PETSC_SUCCESS);
    2737: }

    2739: /* If the file is old, it not only has different path to the coordinates, but   */
    2740: /* does not contain coordinateDMs, so must fall back to the old implementation. */
    2741: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
    2742: {
    2743:   PetscSection coordSection;
    2744:   Vec          coordinates;
    2745:   PetscReal    lengthScale;
    2746:   PetscInt     spatialDim, N, numVertices, vStart, vEnd, v;
    2747:   PetscMPIInt  rank;

    2749:   PetscFunctionBegin;
    2750:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
    2751:   /* Read geometry */
    2752:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
    2753:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates));
    2754:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "vertices"));
    2755:   {
    2756:     /* Force serial load */
    2757:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N));
    2758:     PetscCall(VecSetSizes(coordinates, rank == 0 ? N : 0, N));
    2759:     PetscCall(VecSetBlockSize(coordinates, spatialDim));
    2760:   }
    2761:   PetscCall(VecLoad(coordinates, viewer));
    2762:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    2763:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    2764:   PetscCall(VecScale(coordinates, 1.0 / lengthScale));
    2765:   PetscCall(VecGetLocalSize(coordinates, &numVertices));
    2766:   PetscCall(VecGetBlockSize(coordinates, &spatialDim));
    2767:   numVertices /= spatialDim;
    2768:   /* Create coordinates */
    2769:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    2770:   PetscCheck(numVertices == vEnd - vStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %" PetscInt_FMT " does not match number of vertices %" PetscInt_FMT, numVertices, vEnd - vStart);
    2771:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
    2772:   PetscCall(PetscSectionSetNumFields(coordSection, 1));
    2773:   PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spatialDim));
    2774:   PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd));
    2775:   for (v = vStart; v < vEnd; ++v) {
    2776:     PetscCall(PetscSectionSetDof(coordSection, v, spatialDim));
    2777:     PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spatialDim));
    2778:   }
    2779:   PetscCall(PetscSectionSetUp(coordSection));
    2780:   PetscCall(DMSetCoordinates(dm, coordinates));
    2781:   PetscCall(VecDestroy(&coordinates));
    2782:   PetscFunctionReturn(PETSC_SUCCESS);
    2783: }

    2785: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
    2786: {
    2787:   DMPlexStorageVersion version;
    2788:   DM                   cdm;
    2789:   Vec                  coords;
    2790:   PetscInt             blockSize;
    2791:   PetscReal            lengthScale;
    2792:   PetscSF              lsf;
    2793:   const char          *topologydm_name;
    2794:   char                *coordinatedm_name, *coordinates_name;

    2796:   PetscFunctionBegin;
    2797:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
    2798:   if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
    2799:     PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
    2800:     PetscFunctionReturn(PETSC_SUCCESS);
    2801:   }
    2802:   /* else: since DMPlexStorageVersion 2.0.0 */
    2803:   PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
    2804:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    2805:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    2806:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    2807:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name));
    2808:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name));
    2809:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    2810:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    2811:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    2812:   PetscCall(PetscObjectSetName((PetscObject)cdm, coordinatedm_name));
    2813:   PetscCall(PetscFree(coordinatedm_name));
    2814:   /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
    2815:   PetscCall(DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf));
    2816:   PetscCall(DMCreateLocalVector(cdm, &coords));
    2817:   PetscCall(PetscObjectSetName((PetscObject)coords, coordinates_name));
    2818:   PetscCall(PetscFree(coordinates_name));
    2819:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
    2820:   PetscCall(DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords));
    2821:   PetscCall(PetscViewerPopFormat(viewer));
    2822:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
    2823:   PetscCall(VecScale(coords, 1.0 / lengthScale));
    2824:   PetscCall(DMSetCoordinatesLocal(dm, coords));
    2825:   PetscCall(VecGetBlockSize(coords, &blockSize));
    2826:   PetscCall(DMSetCoordinateDim(dm, blockSize));
    2827:   PetscCall(VecDestroy(&coords));
    2828:   PetscCall(PetscSFDestroy(&lsf));
    2829:   PetscFunctionReturn(PETSC_SUCCESS);
    2830: }

    2832: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
    2833: {
    2834:   DMPlexStorageVersion version;

    2836:   PetscFunctionBegin;
    2837:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
    2838:   PetscCall(PetscInfo(dm, "Loading DM %s storage version %d.%d.%d\n", dm->hdr.name, version->major, version->minor, version->subminor));
    2839:   if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
    2840:     PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL));
    2841:     PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL));
    2842:     PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
    2843:   } else {
    2844:     PetscSF sfXC;

    2846:     /* since DMPlexStorageVersion 2.0.0 */
    2847:     PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC));
    2848:     PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC));
    2849:     PetscCall(DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC));
    2850:     PetscCall(PetscSFDestroy(&sfXC));
    2851:   }
    2852:   PetscFunctionReturn(PETSC_SUCCESS);
    2853: }

    2855: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
    2856: {
    2857:   MPI_Comm  comm;
    2858:   PetscInt  pStart, pEnd, p, m;
    2859:   PetscInt *goffs, *ilocal;
    2860:   PetscBool rootIncludeConstraints, leafIncludeConstraints;

    2862:   PetscFunctionBegin;
    2863:   PetscCall(PetscObjectGetComm((PetscObject)leafSection, &comm));
    2864:   PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
    2865:   PetscCall(PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints));
    2866:   PetscCall(PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints));
    2867:   if (rootIncludeConstraints && leafIncludeConstraints) PetscCall(PetscSectionGetStorageSize(leafSection, &m));
    2868:   else PetscCall(PetscSectionGetConstrainedStorageSize(leafSection, &m));
    2869:   PetscCall(PetscMalloc1(m, &ilocal));
    2870:   PetscCall(PetscMalloc1(m, &goffs));
    2871:   /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
    2872:   /* for the top-level section (not for each field), so one must have   */
    2873:   /* rootSection->pointMajor == PETSC_TRUE.                             */
    2874:   PetscCheck(rootSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
    2875:   /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
    2876:   PetscCheck(leafSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
    2877:   for (p = pStart, m = 0; p < pEnd; ++p) {
    2878:     PetscInt        dof, cdof, i, j, off, goff;
    2879:     const PetscInt *cinds;

    2881:     PetscCall(PetscSectionGetDof(leafSection, p, &dof));
    2882:     if (dof < 0) continue;
    2883:     goff = globalOffsets[p - pStart];
    2884:     PetscCall(PetscSectionGetOffset(leafSection, p, &off));
    2885:     PetscCall(PetscSectionGetConstraintDof(leafSection, p, &cdof));
    2886:     PetscCall(PetscSectionGetConstraintIndices(leafSection, p, &cinds));
    2887:     for (i = 0, j = 0; i < dof; ++i) {
    2888:       PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]);

    2890:       if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {
    2891:         ilocal[m]  = off++;
    2892:         goffs[m++] = goff++;
    2893:       } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
    2894:       else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
    2895:       if (constrained) ++j;
    2896:     }
    2897:   }
    2898:   PetscCall(PetscSFCreate(comm, sectionSF));
    2899:   PetscCall(PetscSFSetFromOptions(*sectionSF));
    2900:   PetscCall(PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs));
    2901:   PetscCall(PetscFree(goffs));
    2902:   PetscFunctionReturn(PETSC_SUCCESS);
    2903: }

    2905: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
    2906: {
    2907:   MPI_Comm     comm;
    2908:   PetscMPIInt  size, rank;
    2909:   const char  *topologydm_name;
    2910:   const char  *sectiondm_name;
    2911:   PetscSection sectionA, sectionB;
    2912:   PetscBool    has;
    2913:   PetscInt     nX, n, i;
    2914:   PetscSF      sfAB;

    2916:   PetscFunctionBegin;
    2917:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    2918:   PetscCallMPI(MPI_Comm_size(comm, &size));
    2919:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    2920:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    2921:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
    2922:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    2923:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    2924:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
    2925:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
    2926:   /* A: on-disk points                        */
    2927:   /* X: list of global point numbers, [0, NX) */
    2928:   /* B: plex points                           */
    2929:   /* Load raw section (sectionA)              */
    2930:   PetscCall(PetscSectionCreate(comm, &sectionA));
    2931:   PetscCall(PetscViewerHDF5HasGroup(viewer, "section", &has));
    2932:   if (has) PetscCall(PetscSectionLoad(sectionA, viewer));
    2933:   else {
    2934:     // TODO If section is missing, create the default affine section with dim dofs on each vertex. Use PetscSplitOwnership() to split vertices.
    2935:     //   How do I know the total number of vertices?
    2936:     PetscInt dim, Nf = 1, Nv, nv = PETSC_DECIDE;

    2938:     PetscCall(DMGetDimension(dm, &dim));
    2939:     PetscCall(DMPlexGetDepthStratumGlobalSize(dm, 0, &Nv));
    2940:     PetscCall(PetscSectionSetNumFields(sectionA, Nf));
    2941:     PetscCall(PetscSectionSetFieldName(sectionA, 0, "Cartesian"));
    2942:     PetscCall(PetscSectionSetFieldComponents(sectionA, 0, dim));
    2943:     for (PetscInt c = 0; c < dim; ++c) {
    2944:       char axis = 'X' + (char)c;

    2946:       PetscCall(PetscSectionSetComponentName(sectionA, 0, c, &axis));
    2947:     }
    2948:     PetscCall(PetscSplitOwnership(comm, &nv, &Nv));
    2949:     PetscCall(PetscSectionSetChart(sectionA, 0, nv));
    2950:     for (PetscInt p = 0; p < nv; ++p) {
    2951:       PetscCall(PetscSectionSetDof(sectionA, p, dim));
    2952:       PetscCall(PetscSectionSetFieldDof(sectionA, p, 0, dim));
    2953:     }
    2954:     PetscCall(PetscSectionSetUp(sectionA));
    2955:   }
    2956:   PetscCall(PetscSectionGetChart(sectionA, NULL, &n));
    2957: /* Create sfAB: A -> B */
    2958: #if defined(PETSC_USE_DEBUG)
    2959:   {
    2960:     PetscInt N, N1;

    2962:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1));
    2963:     PetscCallMPI(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm));
    2964:     PetscCheck(N1 == N, comm, PETSC_ERR_ARG_SIZ, "Mismatching sizes: on-disk order array size (%" PetscInt_FMT ") != number of loaded section points (%" PetscInt_FMT ")", N1, N);
    2965:   }
    2966: #endif
    2967:   {
    2968:     IS              orderIS;
    2969:     const PetscInt *gpoints;
    2970:     PetscSF         sfXA, sfAX;
    2971:     PetscLayout     layout;
    2972:     PetscSFNode    *owners, *buffer;
    2973:     PetscInt        nleaves;
    2974:     PetscInt       *ilocal;
    2975:     PetscSFNode    *iremote;

    2977:     /* Create sfAX: A -> X */
    2978:     PetscCall(ISCreate(comm, &orderIS));
    2979:     PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
    2980:     PetscCall(PetscLayoutSetLocalSize(orderIS->map, n));
    2981:     PetscCall(ISLoad(orderIS, viewer));
    2982:     PetscCall(PetscLayoutCreate(comm, &layout));
    2983:     PetscCall(PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL));
    2984:     PetscCall(PetscLayoutSetLocalSize(layout, nX));
    2985:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
    2986:     PetscCall(PetscLayoutSetUp(layout));
    2987:     PetscCall(PetscSFCreate(comm, &sfXA));
    2988:     PetscCall(ISGetIndices(orderIS, &gpoints));
    2989:     PetscCall(PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints));
    2990:     PetscCall(ISRestoreIndices(orderIS, &gpoints));
    2991:     PetscCall(ISDestroy(&orderIS));
    2992:     PetscCall(PetscLayoutDestroy(&layout));
    2993:     PetscCall(PetscMalloc1(n, &owners));
    2994:     PetscCall(PetscMalloc1(nX, &buffer));
    2995:     for (i = 0; i < n; ++i) {
    2996:       owners[i].rank  = rank;
    2997:       owners[i].index = i;
    2998:     }
    2999:     for (i = 0; i < nX; ++i) {
    3000:       buffer[i].rank  = -1;
    3001:       buffer[i].index = -1;
    3002:     }
    3003:     PetscCall(PetscSFReduceBegin(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC));
    3004:     PetscCall(PetscSFReduceEnd(sfXA, MPIU_SF_NODE, owners, buffer, MPI_MAXLOC));
    3005:     PetscCall(PetscSFDestroy(&sfXA));
    3006:     PetscCall(PetscFree(owners));
    3007:     for (i = 0, nleaves = 0; i < nX; ++i)
    3008:       if (buffer[i].rank >= 0) nleaves++;
    3009:     PetscCall(PetscMalloc1(nleaves, &ilocal));
    3010:     PetscCall(PetscMalloc1(nleaves, &iremote));
    3011:     for (i = 0, nleaves = 0; i < nX; ++i) {
    3012:       if (buffer[i].rank >= 0) {
    3013:         ilocal[nleaves]        = i;
    3014:         iremote[nleaves].rank  = buffer[i].rank;
    3015:         iremote[nleaves].index = buffer[i].index;
    3016:         nleaves++;
    3017:       }
    3018:     }
    3019:     PetscCall(PetscFree(buffer));
    3020:     PetscCall(PetscSFCreate(comm, &sfAX));
    3021:     PetscCall(PetscSFSetFromOptions(sfAX));
    3022:     PetscCall(PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
    3023:     PetscCall(PetscSFCompose(sfAX, sfXB, &sfAB));
    3024:     PetscCall(PetscSFDestroy(&sfAX));
    3025:   }
    3026:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3027:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3028:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3029:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3030:   /* Create plex section (sectionB) */
    3031:   PetscCall(DMGetLocalSection(sectiondm, &sectionB));
    3032:   if (lsf || gsf) {
    3033:     PetscLayout layout;
    3034:     PetscInt    M, m;
    3035:     PetscInt   *offsetsA;
    3036:     PetscBool   includesConstraintsA;

    3038:     PetscCall(PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB));
    3039:     PetscCall(PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA));
    3040:     if (includesConstraintsA) PetscCall(PetscSectionGetStorageSize(sectionA, &m));
    3041:     else PetscCall(PetscSectionGetConstrainedStorageSize(sectionA, &m));
    3042:     PetscCallMPI(MPIU_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm));
    3043:     PetscCall(PetscLayoutCreate(comm, &layout));
    3044:     PetscCall(PetscLayoutSetSize(layout, M));
    3045:     PetscCall(PetscLayoutSetUp(layout));
    3046:     if (lsf) {
    3047:       PetscSF lsfABdata;

    3049:       PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata));
    3050:       *lsf = lsfABdata;
    3051:     }
    3052:     if (gsf) {
    3053:       PetscSection gsectionB, gsectionB1;
    3054:       PetscBool    includesConstraintsB;
    3055:       PetscSF      gsfABdata, pointsf;

    3057:       PetscCall(DMGetGlobalSection(sectiondm, &gsectionB1));
    3058:       PetscCall(PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB));
    3059:       PetscCall(DMGetPointSF(sectiondm, &pointsf));
    3060:       PetscCall(PetscSectionCreateGlobalSection(sectionB, pointsf, PETSC_TRUE, includesConstraintsB, PETSC_TRUE, &gsectionB));
    3061:       PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata));
    3062:       PetscCall(PetscSectionDestroy(&gsectionB));
    3063:       *gsf = gsfABdata;
    3064:     }
    3065:     PetscCall(PetscLayoutDestroy(&layout));
    3066:     PetscCall(PetscFree(offsetsA));
    3067:   } else {
    3068:     PetscCall(PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB));
    3069:   }
    3070:   PetscCall(PetscSFDestroy(&sfAB));
    3071:   PetscCall(PetscSectionDestroy(&sectionA));
    3072:   PetscFunctionReturn(PETSC_SUCCESS);
    3073: }

    3075: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
    3076: {
    3077:   MPI_Comm           comm;
    3078:   const char        *topologydm_name;
    3079:   const char        *sectiondm_name;
    3080:   const char        *vec_name;
    3081:   Vec                vecA;
    3082:   PetscInt           mA, m, bs;
    3083:   const PetscInt    *ilocal;
    3084:   const PetscScalar *src;
    3085:   PetscScalar       *dest;

    3087:   PetscFunctionBegin;
    3088:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    3089:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
    3090:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
    3091:   PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
    3092:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
    3093:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
    3094:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
    3095:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
    3096:   PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
    3097:   PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
    3098:   PetscCall(VecCreate(comm, &vecA));
    3099:   PetscCall(PetscObjectSetName((PetscObject)vecA, vec_name));
    3100:   PetscCall(PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL));
    3101:   /* Check consistency */
    3102:   {
    3103:     PetscSF  pointsf, pointsf1;
    3104:     PetscInt m1, i, j;

    3106:     PetscCall(DMGetPointSF(dm, &pointsf));
    3107:     PetscCall(DMGetPointSF(sectiondm, &pointsf1));
    3108:     PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
    3109: #if defined(PETSC_USE_DEBUG)
    3110:     {
    3111:       PetscInt MA, MA1;

    3113:       PetscCallMPI(MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm));
    3114:       PetscCall(PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1));
    3115:       PetscCheck(MA1 == MA, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%" PetscInt_FMT ") != On-disk vector data size (%" PetscInt_FMT ")", MA, MA1);
    3116:     }
    3117: #endif
    3118:     PetscCall(VecGetLocalSize(vec, &m1));
    3119:     PetscCheck(m1 >= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%" PetscInt_FMT ") < SF leaf size (%" PetscInt_FMT ")", m1, m);
    3120:     for (i = 0; i < m; ++i) {
    3121:       j = ilocal ? ilocal[i] : i;
    3122:       PetscCheck(j >= 0 && j < m1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Leaf's %" PetscInt_FMT "-th index, %" PetscInt_FMT ", not in [%" PetscInt_FMT ", %" PetscInt_FMT ")", i, j, (PetscInt)0, m1);
    3123:     }
    3124:   }
    3125:   PetscCall(VecSetSizes(vecA, mA, PETSC_DECIDE));
    3126:   PetscCall(VecLoad(vecA, viewer));
    3127:   PetscCall(VecGetArrayRead(vecA, &src));
    3128:   PetscCall(VecGetArray(vec, &dest));
    3129:   PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
    3130:   PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
    3131:   PetscCall(VecRestoreArray(vec, &dest));
    3132:   PetscCall(VecRestoreArrayRead(vecA, &src));
    3133:   PetscCall(VecDestroy(&vecA));
    3134:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs));
    3135:   PetscCall(VecSetBlockSize(vec, bs));
    3136:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3137:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3138:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3139:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3140:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3141:   PetscCall(PetscViewerHDF5PopGroup(viewer));
    3142:   PetscFunctionReturn(PETSC_SUCCESS);
    3143: }