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: /* Logging support */
  8: PetscLogEvent DMPLEX_DistributionView, DMPLEX_DistributionLoad;

 10: #if defined(PETSC_HAVE_HDF5)
 11: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);
 12: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer, DMPlexStorageVersion);
 13: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer, const char[], DMPlexStorageVersion);
 14: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer, const char[], DMPlexStorageVersion *);

 16: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

 18: static PetscErrorCode PetscViewerParseVersion_Private(PetscViewer viewer, const char str[], DMPlexStorageVersion *version)
 19: {
 20:   PetscToken           t;
 21:   char                *ts;
 22:   PetscInt             i;
 23:   PetscInt             ti[3];
 24:   DMPlexStorageVersion v;

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

 45: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion v)
 46: {
 47:   PetscContainer cont;

 49:   PetscFunctionBegin;
 50:   PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)viewer), &cont));
 51:   PetscCall(PetscContainerSetPointer(cont, v));
 52:   PetscCall(PetscContainerSetUserDestroy(cont, PetscContainerUserDestroyDefault));
 53:   PetscCall(PetscObjectCompose((PetscObject)viewer, key, (PetscObject)cont));
 54:   PetscCall(PetscContainerDestroy(&cont));
 55:   PetscFunctionReturn(PETSC_SUCCESS);
 56: }

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

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

 69: /*
 70:   Version log:
 71:   1.0.0 legacy version (default if no "dmplex_storage_version" attribute found in file)
 72:   2.0.0 introduce versioning and multiple topologies
 73:   2.1.0 introduce distributions
 74: */
 75: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer viewer, DMPlexStorageVersion version)
 76: {
 77:   PetscBool valid = PETSC_FALSE;

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

126: static inline PetscBool DMPlexStorageVersionGE(DMPlexStorageVersion version, int major, int minor, int subminor)
127: {
128:   return (PetscBool)((version->major == major && version->minor == minor && version->subminor >= subminor) || (version->major == major && version->minor > minor) || (version->major > major));
129: }

131: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion *version)
132: {
133:   const char ATTR_NAME[] = "dmplex_storage_version";
134:   PetscBool  fileHasVersion;
135:   char       optVersion[16], fileVersion[16];

137:   PetscFunctionBegin;
140:   PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, version));
141:   if (*version) PetscFunctionReturn(PETSC_SUCCESS);

143:   PetscCall(PetscStrncpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE, sizeof(fileVersion)));
144:   PetscCall(PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion));
145:   if (fileHasVersion) {
146:     char *tmp;

148:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp));
149:     PetscCall(PetscStrncpy(fileVersion, tmp, sizeof(fileVersion)));
150:     PetscCall(PetscFree(tmp));
151:   }
152:   PetscCall(PetscStrncpy(optVersion, fileVersion, sizeof(optVersion)));
153:   PetscObjectOptionsBegin((PetscObject)viewer);
154:   PetscCall(PetscOptionsString("-dm_plex_view_hdf5_storage_version", "DMPlex HDF5 viewer storage version", NULL, optVersion, optVersion, sizeof(optVersion), NULL));
155:   PetscOptionsEnd();
156:   if (!fileHasVersion) {
157:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, optVersion));
158:   } else {
159:     PetscBool flg;

161:     PetscCall(PetscStrcmp(fileVersion, optVersion, &flg));
162:     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);
163:   }
164:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT));
165:   PetscCall(PetscViewerParseVersion_Private(viewer, optVersion, version));
166:   PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, *version));
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }

170: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion *version)
171: {
172:   const char ATTR_NAME[] = "dmplex_storage_version";
173:   char      *defaultVersion;
174:   char      *versionString;

176:   PetscFunctionBegin;
179:   PetscCall(PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, version));
180:   if (*version) PetscFunctionReturn(PETSC_SUCCESS);

182:   //TODO string HDF5 attribute handling is terrible and should be redesigned
183:   PetscCall(PetscStrallocpy(DMPLEX_STORAGE_VERSION_FIRST, &defaultVersion));
184:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString));
185:   PetscCall(PetscViewerParseVersion_Private(viewer, versionString, version));
186:   PetscCall(PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, *version));
187:   PetscCall(PetscFree(versionString));
188:   PetscCall(PetscFree(defaultVersion));
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

192: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
193: {
194:   PetscFunctionBegin;
195:   if (((PetscObject)dm)->name) {
196:     PetscCall(PetscObjectGetName((PetscObject)dm, name));
197:   } else {
198:     *name = "plex";
199:   }
200:   PetscFunctionReturn(PETSC_SUCCESS);
201: }

203: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
204: {
205:   Vec         stamp;
206:   PetscMPIInt rank;

208:   PetscFunctionBegin;
209:   if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS);
210:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
211:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp));
212:   PetscCall(VecSetBlockSize(stamp, 1));
213:   PetscCall(PetscObjectSetName((PetscObject)stamp, seqname));
214:   if (rank == 0) {
215:     PetscReal timeScale;
216:     PetscBool istime;

218:     PetscCall(PetscStrncmp(seqname, "time", 5, &istime));
219:     if (istime) {
220:       PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale));
221:       value *= timeScale;
222:     }
223:     PetscCall(VecSetValue(stamp, 0, value, INSERT_VALUES));
224:   }
225:   PetscCall(VecAssemblyBegin(stamp));
226:   PetscCall(VecAssemblyEnd(stamp));
227:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/"));
228:   PetscCall(PetscViewerHDF5PushTimestepping(viewer));
229:   PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */
230:   PetscCall(VecView(stamp, viewer));
231:   PetscCall(PetscViewerHDF5PopTimestepping(viewer));
232:   PetscCall(PetscViewerHDF5PopGroup(viewer));
233:   PetscCall(VecDestroy(&stamp));
234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }

237: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
238: {
239:   Vec         stamp;
240:   PetscMPIInt rank;

242:   PetscFunctionBegin;
243:   if (seqnum < 0) PetscFunctionReturn(PETSC_SUCCESS);
244:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank));
245:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp));
246:   PetscCall(VecSetBlockSize(stamp, 1));
247:   PetscCall(PetscObjectSetName((PetscObject)stamp, seqname));
248:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/"));
249:   PetscCall(PetscViewerHDF5PushTimestepping(viewer));
250:   PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum)); /* seqnum < 0 jumps out above */
251:   PetscCall(VecLoad(stamp, viewer));
252:   PetscCall(PetscViewerHDF5PopTimestepping(viewer));
253:   PetscCall(PetscViewerHDF5PopGroup(viewer));
254:   if (rank == 0) {
255:     const PetscScalar *a;
256:     PetscReal          timeScale;
257:     PetscBool          istime;

259:     PetscCall(VecGetArrayRead(stamp, &a));
260:     *value = a[0];
261:     PetscCall(VecRestoreArrayRead(stamp, &a));
262:     PetscCall(PetscStrncmp(seqname, "time", 5, &istime));
263:     if (istime) {
264:       PetscCall(DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale));
265:       *value /= timeScale;
266:     }
267:   }
268:   PetscCall(VecDestroy(&stamp));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
273: {
274:   IS              cutcells = NULL;
275:   const PetscInt *cutc;
276:   PetscInt        cellHeight, vStart, vEnd, cStart, cEnd, c;

278:   PetscFunctionBegin;
279:   if (!cutLabel) PetscFunctionReturn(PETSC_SUCCESS);
280:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
281:   PetscCall(DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd));
282:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
283:   /* Label vertices that should be duplicated */
284:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel));
285:   PetscCall(DMLabelGetStratumIS(cutLabel, 2, &cutcells));
286:   if (cutcells) {
287:     PetscInt n;

289:     PetscCall(ISGetIndices(cutcells, &cutc));
290:     PetscCall(ISGetLocalSize(cutcells, &n));
291:     for (c = 0; c < n; ++c) {
292:       if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
293:         PetscInt *closure = NULL;
294:         PetscInt  closureSize, cl, value;

296:         PetscCall(DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure));
297:         for (cl = 0; cl < closureSize * 2; cl += 2) {
298:           if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
299:             PetscCall(DMLabelGetValue(cutLabel, closure[cl], &value));
300:             if (value == 1) PetscCall(DMLabelSetValue(*cutVertexLabel, closure[cl], 1));
301:           }
302:         }
303:         PetscCall(DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure));
304:       }
305:     }
306:     PetscCall(ISRestoreIndices(cutcells, &cutc));
307:   }
308:   PetscCall(ISDestroy(&cutcells));
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
313: {
314:   DM                      dm;
315:   DM                      dmBC;
316:   PetscSection            section, sectionGlobal;
317:   Vec                     gv;
318:   const char             *name;
319:   PetscViewerVTKFieldType ft;
320:   PetscViewerFormat       format;
321:   PetscInt                seqnum;
322:   PetscReal               seqval;
323:   PetscBool               isseq;

325:   PetscFunctionBegin;
326:   PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
327:   PetscCall(VecGetDM(v, &dm));
328:   PetscCall(DMGetLocalSection(dm, &section));
329:   PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, &seqval));
330:   PetscCall(DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar)seqval, viewer));
331:   if (seqnum >= 0) {
332:     PetscCall(PetscViewerHDF5PushTimestepping(viewer));
333:     PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
334:   }
335:   PetscCall(PetscViewerGetFormat(viewer, &format));
336:   PetscCall(DMGetOutputDM(dm, &dmBC));
337:   PetscCall(DMGetGlobalSection(dmBC, &sectionGlobal));
338:   PetscCall(DMGetGlobalVector(dmBC, &gv));
339:   PetscCall(PetscObjectGetName((PetscObject)v, &name));
340:   PetscCall(PetscObjectSetName((PetscObject)gv, name));
341:   PetscCall(DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv));
342:   PetscCall(DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv));
343:   PetscCall(PetscObjectTypeCompare((PetscObject)gv, VECSEQ, &isseq));
344:   if (isseq) PetscCall(VecView_Seq(gv, viewer));
345:   else PetscCall(VecView_MPI(gv, viewer));
346:   if (format == PETSC_VIEWER_HDF5_VIZ) {
347:     /* Output visualization representation */
348:     PetscInt numFields, f;
349:     DMLabel  cutLabel, cutVertexLabel = NULL;

351:     PetscCall(PetscSectionGetNumFields(section, &numFields));
352:     PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
353:     for (f = 0; f < numFields; ++f) {
354:       Vec         subv;
355:       IS          is;
356:       const char *fname, *fgroup, *componentName;
357:       char        subname[PETSC_MAX_PATH_LEN];
358:       PetscInt    pStart, pEnd, Nc, c;

360:       PetscCall(DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft));
361:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
362:       PetscCall(PetscSectionGetFieldName(section, f, &fname));
363:       if (!fname || ft == PETSC_VTK_INVALID) continue;
364:       PetscCall(PetscViewerHDF5PushGroup(viewer, fgroup));
365:       if (cutLabel) {
366:         const PetscScalar *ga;
367:         PetscScalar       *suba;
368:         PetscInt           gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;

370:         PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
371:         PetscCall(PetscSectionGetFieldComponents(section, f, &Nc));
372:         for (p = pStart; p < pEnd; ++p) {
373:           PetscInt gdof, fdof = 0, val;

375:           PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
376:           if (gdof > 0) PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
377:           subSize += fdof;
378:           PetscCall(DMLabelGetValue(cutVertexLabel, p, &val));
379:           if (val == 1) extSize += fdof;
380:         }
381:         PetscCall(VecCreate(PetscObjectComm((PetscObject)gv), &subv));
382:         PetscCall(VecSetSizes(subv, subSize + extSize, PETSC_DETERMINE));
383:         PetscCall(VecSetBlockSize(subv, Nc));
384:         PetscCall(VecSetType(subv, VECSTANDARD));
385:         PetscCall(VecGetOwnershipRange(gv, &gstart, NULL));
386:         PetscCall(VecGetArrayRead(gv, &ga));
387:         PetscCall(VecGetArray(subv, &suba));
388:         for (p = pStart; p < pEnd; ++p) {
389:           PetscInt gdof, goff, val;

391:           PetscCall(PetscSectionGetDof(sectionGlobal, p, &gdof));
392:           if (gdof > 0) {
393:             PetscInt fdof, fc, f2, poff = 0;

395:             PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff));
396:             /* Can get rid of this loop by storing field information in the global section */
397:             for (f2 = 0; f2 < f; ++f2) {
398:               PetscCall(PetscSectionGetFieldDof(section, p, f2, &fdof));
399:               poff += fdof;
400:             }
401:             PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
402:             for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff + poff + fc - gstart];
403:             PetscCall(DMLabelGetValue(cutVertexLabel, p, &val));
404:             if (val == 1) {
405:               for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize + newOff] = ga[goff + poff + fc - gstart];
406:             }
407:           }
408:         }
409:         PetscCall(VecRestoreArrayRead(gv, &ga));
410:         PetscCall(VecRestoreArray(subv, &suba));
411:         PetscCall(DMLabelDestroy(&cutVertexLabel));
412:       } else {
413:         PetscCall(PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv));
414:       }
415:       PetscCall(PetscStrncpy(subname, name, sizeof(subname)));
416:       PetscCall(PetscStrlcat(subname, "_", sizeof(subname)));
417:       PetscCall(PetscStrlcat(subname, fname, sizeof(subname)));
418:       PetscCall(PetscObjectSetName((PetscObject)subv, subname));
419:       if (isseq) PetscCall(VecView_Seq(subv, viewer));
420:       else PetscCall(VecView_MPI(subv, viewer));
421:       if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
422:         PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "vector"));
423:       } else {
424:         PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "scalar"));
425:       }

427:       /* Output the component names in the field if available */
428:       PetscCall(PetscSectionGetFieldComponents(section, f, &Nc));
429:       for (c = 0; c < Nc; ++c) {
430:         char componentNameLabel[PETSC_MAX_PATH_LEN];
431:         PetscCall(PetscSectionGetComponentName(section, f, c, &componentName));
432:         PetscCall(PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%" PetscInt_FMT, c));
433:         PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, componentNameLabel, PETSC_STRING, componentName));
434:       }

436:       if (cutLabel) PetscCall(VecDestroy(&subv));
437:       else PetscCall(PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv));
438:       PetscCall(PetscViewerHDF5PopGroup(viewer));
439:     }
440:   }
441:   if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
442:   PetscCall(DMRestoreGlobalVector(dmBC, &gv));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
447: {
448:   DM          dm;
449:   Vec         locv;
450:   PetscObject isZero;
451:   const char *name;
452:   PetscReal   time;

454:   PetscFunctionBegin;
455:   PetscCall(VecGetDM(v, &dm));
456:   PetscCall(DMGetLocalVector(dm, &locv));
457:   PetscCall(PetscObjectGetName((PetscObject)v, &name));
458:   PetscCall(PetscObjectSetName((PetscObject)locv, name));
459:   PetscCall(PetscObjectQuery((PetscObject)v, "__Vec_bc_zero__", &isZero));
460:   PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", isZero));
461:   PetscCall(DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv));
462:   PetscCall(DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv));
463:   PetscCall(DMGetOutputSequenceNumber(dm, NULL, &time));
464:   PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL));
465:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
466:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
467:   PetscCall(VecView_Plex_Local_HDF5_Internal(locv, viewer));
468:   PetscCall(PetscViewerPopFormat(viewer));
469:   PetscCall(PetscViewerHDF5PopGroup(viewer));
470:   PetscCall(PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", NULL));
471:   PetscCall(DMRestoreLocalVector(dm, &locv));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
476: {
477:   PetscBool isseq;

479:   PetscFunctionBegin;
480:   PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq));
481:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
482:   if (isseq) PetscCall(VecView_Seq(v, viewer));
483:   else PetscCall(VecView_MPI(v, viewer));
484:   PetscCall(PetscViewerHDF5PopGroup(viewer));
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

488: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
489: {
490:   DM          dm;
491:   Vec         locv;
492:   const char *name;
493:   PetscInt    seqnum;

495:   PetscFunctionBegin;
496:   PetscCall(VecGetDM(v, &dm));
497:   PetscCall(DMGetLocalVector(dm, &locv));
498:   PetscCall(PetscObjectGetName((PetscObject)v, &name));
499:   PetscCall(PetscObjectSetName((PetscObject)locv, name));
500:   PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL));
501:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
502:   if (seqnum >= 0) {
503:     PetscCall(PetscViewerHDF5PushTimestepping(viewer));
504:     PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
505:   }
506:   PetscCall(VecLoad_Plex_Local(locv, viewer));
507:   if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
508:   PetscCall(PetscViewerHDF5PopGroup(viewer));
509:   PetscCall(DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v));
510:   PetscCall(DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v));
511:   PetscCall(DMRestoreLocalVector(dm, &locv));
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

515: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
516: {
517:   DM       dm;
518:   PetscInt seqnum;

520:   PetscFunctionBegin;
521:   PetscCall(VecGetDM(v, &dm));
522:   PetscCall(DMGetOutputSequenceNumber(dm, &seqnum, NULL));
523:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/fields"));
524:   if (seqnum >= 0) {
525:     PetscCall(PetscViewerHDF5PushTimestepping(viewer));
526:     PetscCall(PetscViewerHDF5SetTimestep(viewer, seqnum));
527:   }
528:   PetscCall(VecLoad_Default(v, viewer));
529:   if (seqnum >= 0) PetscCall(PetscViewerHDF5PopTimestepping(viewer));
530:   PetscCall(PetscViewerHDF5PopGroup(viewer));
531:   PetscFunctionReturn(PETSC_SUCCESS);
532: }

534: static PetscErrorCode DMPlexDistributionView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
535: {
536:   MPI_Comm           comm;
537:   PetscMPIInt        size, rank;
538:   PetscInt           size_petsc_int;
539:   const char        *topologydm_name, *distribution_name;
540:   const PetscInt    *gpoint;
541:   PetscInt           pStart, pEnd, p;
542:   PetscSF            pointSF;
543:   PetscInt           nroots, nleaves;
544:   const PetscInt    *ilocal;
545:   const PetscSFNode *iremote;
546:   IS                 chartSizesIS, ownersIS, gpointsIS;
547:   PetscInt          *chartSize, *owners, *gpoints;

549:   PetscFunctionBegin;
550:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
551:   PetscCallMPI(MPI_Comm_size(comm, &size));
552:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
553:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
554:   PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
555:   if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
556:   PetscCall(PetscLogEventBegin(DMPLEX_DistributionView, viewer, 0, 0, 0));
557:   size_petsc_int = (PetscInt)size;
558:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size_petsc_int));
559:   PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
560:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
561:   PetscCall(PetscMalloc1(1, &chartSize));
562:   *chartSize = pEnd - pStart;
563:   PetscCall(PetscMalloc1(*chartSize, &owners));
564:   PetscCall(PetscMalloc1(*chartSize, &gpoints));
565:   PetscCall(DMGetPointSF(dm, &pointSF));
566:   PetscCall(PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote));
567:   for (p = pStart; p < pEnd; ++p) {
568:     PetscInt gp = gpoint[p - pStart];

570:     owners[p - pStart]  = rank;
571:     gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp);
572:   }
573:   for (p = 0; p < nleaves; ++p) {
574:     PetscInt ilocalp = (ilocal ? ilocal[p] : p);

576:     owners[ilocalp] = iremote[p].rank;
577:   }
578:   PetscCall(ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS));
579:   PetscCall(ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS));
580:   PetscCall(ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS));
581:   PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
582:   PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
583:   PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
584:   PetscCall(ISView(chartSizesIS, viewer));
585:   PetscCall(ISView(ownersIS, viewer));
586:   PetscCall(ISView(gpointsIS, viewer));
587:   PetscCall(ISDestroy(&chartSizesIS));
588:   PetscCall(ISDestroy(&ownersIS));
589:   PetscCall(ISDestroy(&gpointsIS));
590:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
591:   PetscCall(PetscLogEventEnd(DMPLEX_DistributionView, viewer, 0, 0, 0));
592:   PetscFunctionReturn(PETSC_SUCCESS);
593: }

595: 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[])
596: {
597:   IS              coneSizesIS, conesIS, orientationsIS;
598:   PetscInt       *coneSizes, *cones, *orientations;
599:   const PetscInt *gpoint;
600:   PetscInt        nPoints = 0, conesSize = 0;
601:   PetscInt        p, c, s;
602:   MPI_Comm        comm;

604:   PetscFunctionBegin;
605:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
606:   PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
607:   for (p = pStart; p < pEnd; ++p) {
608:     if (gpoint[p] >= 0) {
609:       PetscInt coneSize;

611:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
612:       nPoints += 1;
613:       conesSize += coneSize;
614:     }
615:   }
616:   PetscCall(PetscMalloc1(nPoints, &coneSizes));
617:   PetscCall(PetscMalloc1(conesSize, &cones));
618:   PetscCall(PetscMalloc1(conesSize, &orientations));
619:   for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
620:     if (gpoint[p] >= 0) {
621:       const PetscInt *cone, *ornt;
622:       PetscInt        coneSize, cp;

624:       PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
625:       PetscCall(DMPlexGetCone(dm, p, &cone));
626:       PetscCall(DMPlexGetConeOrientation(dm, p, &ornt));
627:       coneSizes[s] = coneSize;
628:       for (cp = 0; cp < coneSize; ++cp, ++c) {
629:         cones[c]        = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]];
630:         orientations[c] = ornt[cp];
631:       }
632:       ++s;
633:     }
634:   }
635:   PetscCheck(s == nPoints, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %" PetscInt_FMT " != %" PetscInt_FMT, s, nPoints);
636:   PetscCheck(c == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %" PetscInt_FMT " != %" PetscInt_FMT, c, conesSize);
637:   PetscCall(ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS));
638:   PetscCall(ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS));
639:   PetscCall(ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS));
640:   PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
641:   PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
642:   PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
643:   PetscCall(ISView(coneSizesIS, viewer));
644:   PetscCall(ISView(conesIS, viewer));
645:   PetscCall(ISView(orientationsIS, viewer));
646:   PetscCall(ISDestroy(&coneSizesIS));
647:   PetscCall(ISDestroy(&conesIS));
648:   PetscCall(ISDestroy(&orientationsIS));
649:   if (pointsName) {
650:     IS        pointsIS;
651:     PetscInt *points;

653:     PetscCall(PetscMalloc1(nPoints, &points));
654:     for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
655:       if (gpoint[p] >= 0) {
656:         points[s] = gpoint[p];
657:         ++s;
658:       }
659:     }
660:     PetscCall(ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS));
661:     PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
662:     PetscCall(ISView(pointsIS, viewer));
663:     PetscCall(ISDestroy(&pointsIS));
664:   }
665:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
666:   PetscFunctionReturn(PETSC_SUCCESS);
667: }

669: static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
670: {
671:   const char *pointsName, *coneSizesName, *conesName, *orientationsName;
672:   PetscInt    pStart, pEnd, dim;

674:   PetscFunctionBegin;
675:   pointsName       = "order";
676:   coneSizesName    = "cones";
677:   conesName        = "cells";
678:   orientationsName = "orientation";
679:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
680:   PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName));
681:   PetscCall(DMGetDimension(dm, &dim));
682:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim));
683:   PetscFunctionReturn(PETSC_SUCCESS);
684: }

686: //TODO get this numbering right away without needing this function
687: /* Renumber global point numbers so that they are 0-based per stratum */
688: static PetscErrorCode RenumberGlobalPointNumbersPerStratum_Private(DM dm, IS globalPointNumbers, IS *newGlobalPointNumbers, IS *strataPermutation)
689: {
690:   PetscInt        d, depth, p, n;
691:   PetscInt       *offsets;
692:   const PetscInt *gpn;
693:   PetscInt       *ngpn;
694:   MPI_Comm        comm;

696:   PetscFunctionBegin;
697:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
698:   PetscCall(ISGetLocalSize(globalPointNumbers, &n));
699:   PetscCall(ISGetIndices(globalPointNumbers, &gpn));
700:   PetscCall(PetscMalloc1(n, &ngpn));
701:   PetscCall(DMPlexGetDepth(dm, &depth));
702:   PetscCall(PetscMalloc1(depth + 1, &offsets));
703:   for (d = 0; d <= depth; d++) {
704:     PetscInt pStart, pEnd;

706:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
707:     offsets[d] = PETSC_MAX_INT;
708:     for (p = pStart; p < pEnd; p++) {
709:       if (gpn[p] >= 0 && gpn[p] < offsets[d]) offsets[d] = gpn[p];
710:     }
711:   }
712:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, offsets, depth + 1, MPIU_INT, MPI_MIN, comm));
713:   for (d = 0; d <= depth; d++) {
714:     PetscInt pStart, pEnd;

716:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
717:     for (p = pStart; p < pEnd; p++) ngpn[p] = gpn[p] - PetscSign(gpn[p]) * offsets[d];
718:   }
719:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpn));
720:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)globalPointNumbers), n, ngpn, PETSC_OWN_POINTER, newGlobalPointNumbers));
721:   {
722:     PetscInt *perm;

724:     PetscCall(PetscMalloc1(depth + 1, &perm));
725:     for (d = 0; d <= depth; d++) perm[d] = d;
726:     PetscCall(PetscSortIntWithPermutation(depth + 1, offsets, perm));
727:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, depth + 1, perm, PETSC_OWN_POINTER, strataPermutation));
728:   }
729:   PetscCall(PetscFree(offsets));
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

733: static PetscErrorCode DMPlexTopologyView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
734: {
735:   IS          globalPointNumbers0, strataPermutation;
736:   const char *coneSizesName, *conesName, *orientationsName;
737:   PetscInt    depth, d;
738:   MPI_Comm    comm;

740:   PetscFunctionBegin;
741:   coneSizesName    = "cone_sizes";
742:   conesName        = "cones";
743:   orientationsName = "orientations";
744:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
745:   PetscCall(DMPlexGetDepth(dm, &depth));
746:   {
747:     PetscInt dim;

749:     PetscCall(DMGetDimension(dm, &dim));
750:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "cell_dim", PETSC_INT, &dim));
751:     PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "depth", PETSC_INT, &depth));
752:   }

754:   PetscCall(RenumberGlobalPointNumbersPerStratum_Private(dm, globalPointNumbers, &globalPointNumbers0, &strataPermutation));
755:   /* TODO dirty trick to save serial IS using the same parallel viewer */
756:   {
757:     IS              spOnComm;
758:     PetscInt        n   = 0, N;
759:     const PetscInt *idx = NULL;
760:     const PetscInt *old;
761:     PetscMPIInt     rank;

763:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
764:     PetscCall(ISGetLocalSize(strataPermutation, &N));
765:     PetscCall(ISGetIndices(strataPermutation, &old));
766:     if (!rank) {
767:       n   = N;
768:       idx = old;
769:     }
770:     PetscCall(ISCreateGeneral(comm, n, idx, PETSC_COPY_VALUES, &spOnComm));
771:     PetscCall(ISRestoreIndices(strataPermutation, &old));
772:     PetscCall(ISDestroy(&strataPermutation));
773:     strataPermutation = spOnComm;
774:   }
775:   PetscCall(PetscObjectSetName((PetscObject)strataPermutation, "permutation"));
776:   PetscCall(ISView(strataPermutation, viewer));
777:   PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
778:   for (d = 0; d <= depth; d++) {
779:     PetscInt pStart, pEnd;
780:     char     group[128];

782:     PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, d));
783:     PetscCall(PetscViewerHDF5PushGroup(viewer, group));
784:     PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
785:     PetscCall(DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers0, viewer, pStart, pEnd, NULL, coneSizesName, conesName, orientationsName));
786:     PetscCall(PetscViewerHDF5PopGroup(viewer));
787:   }
788:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */
789:   PetscCall(ISDestroy(&globalPointNumbers0));
790:   PetscCall(ISDestroy(&strataPermutation));
791:   PetscFunctionReturn(PETSC_SUCCESS);
792: }

794: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
795: {
796:   DMPlexStorageVersion version;
797:   const char          *topologydm_name;
798:   char                 group[PETSC_MAX_PATH_LEN];

800:   PetscFunctionBegin;
801:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
802:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
803:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
804:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
805:   } else {
806:     PetscCall(PetscStrncpy(group, "/", sizeof(group)));
807:   }
808:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));

810:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
811:   if (version->major < 3) {
812:     PetscCall(DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer));
813:   } else {
814:     /* since DMPlexStorageVersion 3.0.0 */
815:     PetscCall(DMPlexTopologyView_HDF5_Private(dm, globalPointNumbers, viewer));
816:   }
817:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */

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

822:     PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
823:     PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
824:     PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
825:     PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
826:     PetscCall(DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer));
827:     PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
828:     PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
829:   }

831:   PetscCall(PetscViewerHDF5PopGroup(viewer));
832:   PetscFunctionReturn(PETSC_SUCCESS);
833: }

835: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
836: {
837:   PetscSF         sfPoint;
838:   DMLabel         cutLabel, cutVertexLabel         = NULL;
839:   IS              globalVertexNumbers, cutvertices = NULL;
840:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
841:   PetscInt       *vertices;
842:   PetscInt        conesSize = 0;
843:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;

845:   PetscFunctionBegin;
846:   *numCorners = 0;
847:   PetscCall(DMGetDimension(dm, &dim));
848:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
849:   PetscCall(ISGetIndices(globalCellNumbers, &gcell));

851:   for (cell = cStart; cell < cEnd; ++cell) {
852:     PetscInt *closure = NULL;
853:     PetscInt  closureSize, v, Nc = 0;

855:     if (gcell[cell] < 0) continue;
856:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
857:     for (v = 0; v < closureSize * 2; v += 2) {
858:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
859:     }
860:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
861:     conesSize += Nc;
862:     if (!numCornersLocal) numCornersLocal = Nc;
863:     else if (numCornersLocal != Nc) numCornersLocal = 1;
864:   }
865:   PetscCall(MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
866:   PetscCheck(!numCornersLocal || !(numCornersLocal != *numCorners || *numCorners == 0), PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
867:   /* Handle periodic cuts by identifying vertices which should be duplicated */
868:   PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
869:   PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
870:   if (cutVertexLabel) PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices));
871:   if (cutvertices) {
872:     PetscCall(ISGetIndices(cutvertices, &cutverts));
873:     PetscCall(ISGetLocalSize(cutvertices, &vExtra));
874:   }
875:   PetscCall(DMGetPointSF(dm, &sfPoint));
876:   if (cutLabel) {
877:     const PetscInt    *ilocal;
878:     const PetscSFNode *iremote;
879:     PetscInt           nroots, nleaves;

881:     PetscCall(PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote));
882:     if (nleaves < 0) {
883:       PetscCall(PetscObjectReference((PetscObject)sfPoint));
884:     } else {
885:       PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint));
886:       PetscCall(PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES));
887:     }
888:   } else {
889:     PetscCall(PetscObjectReference((PetscObject)sfPoint));
890:   }
891:   /* Number all vertices */
892:   PetscCall(DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers));
893:   PetscCall(PetscSFDestroy(&sfPoint));
894:   /* Create cones */
895:   PetscCall(ISGetIndices(globalVertexNumbers, &gvertex));
896:   PetscCall(PetscMalloc1(conesSize, &vertices));
897:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
898:     PetscInt *closure = NULL;
899:     PetscInt  closureSize, Nc = 0, p, value = -1;
900:     PetscBool replace;

902:     if (gcell[cell] < 0) continue;
903:     if (cutLabel) PetscCall(DMLabelGetValue(cutLabel, cell, &value));
904:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
905:     PetscCall(DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
906:     for (p = 0; p < closureSize * 2; p += 2) {
907:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p];
908:     }
909:     PetscCall(DMPlexReorderCell(dm, cell, closure));
910:     for (p = 0; p < Nc; ++p) {
911:       PetscInt nv, gv = gvertex[closure[p] - vStart];

913:       if (replace) {
914:         PetscCall(PetscFindInt(closure[p], vExtra, cutverts, &nv));
915:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
916:       }
917:       vertices[v++] = gv < 0 ? -(gv + 1) : gv;
918:     }
919:     PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure));
920:   }
921:   PetscCall(ISRestoreIndices(globalVertexNumbers, &gvertex));
922:   PetscCall(ISDestroy(&globalVertexNumbers));
923:   PetscCall(ISRestoreIndices(globalCellNumbers, &gcell));
924:   if (cutvertices) PetscCall(ISRestoreIndices(cutvertices, &cutverts));
925:   PetscCall(ISDestroy(&cutvertices));
926:   PetscCall(DMLabelDestroy(&cutVertexLabel));
927:   PetscCheck(v == conesSize, PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %" PetscInt_FMT " != %" PetscInt_FMT, v, conesSize);
928:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS));
929:   PetscCall(PetscLayoutSetBlockSize((*cellIS)->map, *numCorners));
930:   PetscCall(PetscObjectSetName((PetscObject)*cellIS, "cells"));
931:   PetscFunctionReturn(PETSC_SUCCESS);
932: }

934: static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer)
935: {
936:   DM       cdm;
937:   DMLabel  depthLabel, ctLabel;
938:   IS       cellIS;
939:   PetscInt dim, depth, cellHeight, c, n = 0;

941:   PetscFunctionBegin;
942:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
943:   PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));

945:   PetscCall(PetscViewerHDF5PopGroup(viewer));
946:   PetscCall(DMGetDimension(dm, &dim));
947:   PetscCall(DMPlexGetDepth(dm, &depth));
948:   PetscCall(DMGetCoordinateDM(dm, &cdm));
949:   PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
950:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
951:   PetscCall(DMPlexGetCellTypeLabel(dm, &ctLabel));
952:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
953:     const DMPolytopeType ict = (DMPolytopeType)c;
954:     PetscInt             pStart, pEnd, dep, numCorners;
955:     PetscBool            output = PETSC_FALSE, doOutput;

957:     if (ict == DM_POLYTOPE_FV_GHOST) continue;
958:     PetscCall(DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd));
959:     if (pStart >= 0) {
960:       PetscCall(DMLabelGetValue(depthLabel, pStart, &dep));
961:       if (dep == depth - cellHeight) output = PETSC_TRUE;
962:     }
963:     PetscCall(MPIU_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm)));
964:     if (!doOutput) continue;
965:     PetscCall(CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS));
966:     if (!n) {
967:       PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/topology"));
968:     } else {
969:       char group[PETSC_MAX_PATH_LEN];

971:       PetscCall(PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n));
972:       PetscCall(PetscViewerHDF5PushGroup(viewer, group));
973:     }
974:     PetscCall(ISView(cellIS, viewer));
975:     PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners));
976:     PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim));
977:     PetscCall(ISDestroy(&cellIS));
978:     PetscCall(PetscViewerHDF5PopGroup(viewer));
979:     ++n;
980:   }
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }

984: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
985: {
986:   DM        cdm;
987:   Vec       coordinates, newcoords;
988:   PetscReal lengthScale;
989:   PetscInt  m, M, bs;

991:   PetscFunctionBegin;
992:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
993:   PetscCall(DMGetCoordinateDM(dm, &cdm));
994:   PetscCall(DMGetCoordinates(dm, &coordinates));
995:   PetscCall(VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords));
996:   PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
997:   PetscCall(VecGetSize(coordinates, &M));
998:   PetscCall(VecGetLocalSize(coordinates, &m));
999:   PetscCall(VecSetSizes(newcoords, m, M));
1000:   PetscCall(VecGetBlockSize(coordinates, &bs));
1001:   PetscCall(VecSetBlockSize(newcoords, bs));
1002:   PetscCall(VecSetType(newcoords, VECSTANDARD));
1003:   PetscCall(VecCopy(coordinates, newcoords));
1004:   PetscCall(VecScale(newcoords, lengthScale));
1005:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
1006:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
1007:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
1008:   PetscCall(VecView(newcoords, viewer));
1009:   PetscCall(PetscViewerPopFormat(viewer));
1010:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1011:   PetscCall(VecDestroy(&newcoords));
1012:   PetscFunctionReturn(PETSC_SUCCESS);
1013: }

1015: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
1016: {
1017:   DM          cdm;
1018:   Vec         coords, newcoords;
1019:   PetscInt    m, M, bs;
1020:   PetscReal   lengthScale;
1021:   const char *topologydm_name, *coordinatedm_name, *coordinates_name;

1023:   PetscFunctionBegin;
1024:   {
1025:     PetscViewerFormat    format;
1026:     DMPlexStorageVersion version;

1028:     PetscCall(PetscViewerGetFormat(viewer, &format));
1029:     PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1030:     if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) {
1031:       PetscCall(DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer));
1032:       PetscFunctionReturn(PETSC_SUCCESS);
1033:     }
1034:   }
1035:   /* since 2.0.0 */
1036:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1037:   PetscCall(DMGetCoordinates(dm, &coords));
1038:   PetscCall(PetscObjectGetName((PetscObject)cdm, &coordinatedm_name));
1039:   PetscCall(PetscObjectGetName((PetscObject)coords, &coordinates_name));
1040:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1041:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1042:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1043:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name));
1044:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name));
1045:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1046:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1047:   PetscCall(DMPlexSectionView(dm, viewer, cdm));
1048:   PetscCall(VecCreate(PetscObjectComm((PetscObject)coords), &newcoords));
1049:   PetscCall(PetscObjectSetName((PetscObject)newcoords, coordinates_name));
1050:   PetscCall(VecGetSize(coords, &M));
1051:   PetscCall(VecGetLocalSize(coords, &m));
1052:   PetscCall(VecSetSizes(newcoords, m, M));
1053:   PetscCall(VecGetBlockSize(coords, &bs));
1054:   PetscCall(VecSetBlockSize(newcoords, bs));
1055:   PetscCall(VecSetType(newcoords, VECSTANDARD));
1056:   PetscCall(VecCopy(coords, newcoords));
1057:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1058:   PetscCall(VecScale(newcoords, lengthScale));
1059:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
1060:   PetscCall(DMPlexGlobalVectorView(dm, viewer, cdm, newcoords));
1061:   PetscCall(PetscViewerPopFormat(viewer));
1062:   PetscCall(VecDestroy(&newcoords));
1063:   PetscFunctionReturn(PETSC_SUCCESS);
1064: }

1066: static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer)
1067: {
1068:   DM               cdm;
1069:   Vec              coordinatesLocal, newcoords;
1070:   PetscSection     cSection, cGlobalSection;
1071:   PetscScalar     *coords, *ncoords;
1072:   DMLabel          cutLabel, cutVertexLabel = NULL;
1073:   const PetscReal *L;
1074:   PetscReal        lengthScale;
1075:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
1076:   PetscBool        localized, embedded;

1078:   PetscFunctionBegin;
1079:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1080:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
1081:   PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
1082:   PetscCall(VecGetBlockSize(coordinatesLocal, &bs));
1083:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
1084:   if (localized == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1085:   PetscCall(DMGetPeriodicity(dm, NULL, NULL, &L));
1086:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1087:   PetscCall(DMGetLocalSection(cdm, &cSection));
1088:   PetscCall(DMGetGlobalSection(cdm, &cGlobalSection));
1089:   PetscCall(DMGetLabel(dm, "periodic_cut", &cutLabel));
1090:   N = 0;

1092:   PetscCall(DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel));
1093:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &newcoords));
1094:   PetscCall(PetscSectionGetDof(cSection, vStart, &dof));
1095:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof));
1096:   embedded = (PetscBool)(L && dof == 2 && !cutLabel);
1097:   if (cutVertexLabel) {
1098:     PetscCall(DMLabelGetStratumSize(cutVertexLabel, 1, &v));
1099:     N += dof * v;
1100:   }
1101:   for (v = vStart; v < vEnd; ++v) {
1102:     PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
1103:     if (dof < 0) continue;
1104:     if (embedded) N += dof + 1;
1105:     else N += dof;
1106:   }
1107:   if (embedded) PetscCall(VecSetBlockSize(newcoords, bs + 1));
1108:   else PetscCall(VecSetBlockSize(newcoords, bs));
1109:   PetscCall(VecSetSizes(newcoords, N, PETSC_DETERMINE));
1110:   PetscCall(VecSetType(newcoords, VECSTANDARD));
1111:   PetscCall(VecGetArray(coordinatesLocal, &coords));
1112:   PetscCall(VecGetArray(newcoords, &ncoords));
1113:   coordSize = 0;
1114:   for (v = vStart; v < vEnd; ++v) {
1115:     PetscCall(PetscSectionGetDof(cGlobalSection, v, &dof));
1116:     PetscCall(PetscSectionGetOffset(cSection, v, &off));
1117:     if (dof < 0) continue;
1118:     if (embedded) {
1119:       if (L && (L[0] > 0.0) && (L[1] > 0.0)) {
1120:         PetscReal theta, phi, r, R;
1121:         /* XY-periodic */
1122:         /* Suppose its an y-z circle, then
1123:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
1124:            and the circle in that plane is
1125:              \hat r cos(phi) + \hat x sin(phi) */
1126:         theta                = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1];
1127:         phi                  = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0];
1128:         r                    = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]);
1129:         R                    = L[1] / (2.0 * PETSC_PI);
1130:         ncoords[coordSize++] = PetscSinReal(phi) * r;
1131:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
1132:         ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
1133:       } else if (L && (L[0] > 0.0)) {
1134:         /* X-periodic */
1135:         ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
1136:         ncoords[coordSize++] = coords[off + 1];
1137:         ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
1138:       } else if (L && (L[1] > 0.0)) {
1139:         /* Y-periodic */
1140:         ncoords[coordSize++] = coords[off + 0];
1141:         ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
1142:         ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
1143:   #if 0
1144:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
1145:         PetscReal phi, r, R;
1146:         /* Mobius strip */
1147:         /* Suppose its an x-z circle, then
1148:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
1149:            and in that plane we rotate by pi as we go around the circle
1150:              \hat r cos(phi/2) + \hat y sin(phi/2) */
1151:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
1152:         R     = L[0];
1153:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
1154:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
1155:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
1156:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
1157:   #endif
1158:       } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
1159:     } else {
1160:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d];
1161:     }
1162:   }
1163:   if (cutVertexLabel) {
1164:     IS              vertices;
1165:     const PetscInt *verts;
1166:     PetscInt        n;

1168:     PetscCall(DMLabelGetStratumIS(cutVertexLabel, 1, &vertices));
1169:     if (vertices) {
1170:       PetscCall(ISGetIndices(vertices, &verts));
1171:       PetscCall(ISGetLocalSize(vertices, &n));
1172:       for (v = 0; v < n; ++v) {
1173:         PetscCall(PetscSectionGetDof(cSection, verts[v], &dof));
1174:         PetscCall(PetscSectionGetOffset(cSection, verts[v], &off));
1175:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0);
1176:       }
1177:       PetscCall(ISRestoreIndices(vertices, &verts));
1178:       PetscCall(ISDestroy(&vertices));
1179:     }
1180:   }
1181:   PetscCheck(coordSize == N, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %" PetscInt_FMT " != %" PetscInt_FMT, coordSize, N);
1182:   PetscCall(DMLabelDestroy(&cutVertexLabel));
1183:   PetscCall(VecRestoreArray(coordinatesLocal, &coords));
1184:   PetscCall(VecRestoreArray(newcoords, &ncoords));
1185:   PetscCall(PetscObjectSetName((PetscObject)newcoords, "vertices"));
1186:   PetscCall(VecScale(newcoords, lengthScale));
1187:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz"));
1188:   PetscCall(PetscViewerHDF5WriteGroup(viewer, NULL));
1189:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1190:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/viz/geometry"));
1191:   PetscCall(VecView(newcoords, viewer));
1192:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1193:   PetscCall(VecDestroy(&newcoords));
1194:   PetscFunctionReturn(PETSC_SUCCESS);
1195: }

1197: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
1198: {
1199:   const char          *topologydm_name;
1200:   const PetscInt      *gpoint;
1201:   PetscInt             numLabels, l;
1202:   DMPlexStorageVersion version;
1203:   char                 group[PETSC_MAX_PATH_LEN];

1205:   PetscFunctionBegin;
1206:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version));
1207:   PetscCall(ISGetIndices(globalPointNumbers, &gpoint));
1208:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1209:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1210:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
1211:   } else {
1212:     PetscCall(PetscStrncpy(group, "/labels", sizeof(group)));
1213:   }
1214:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1215:   PetscCall(DMGetNumLabels(dm, &numLabels));
1216:   for (l = 0; l < numLabels; ++l) {
1217:     DMLabel         label;
1218:     const char     *name;
1219:     IS              valueIS, pvalueIS, globalValueIS;
1220:     const PetscInt *values;
1221:     PetscInt        numValues, v;
1222:     PetscBool       isDepth, output;

1224:     PetscCall(DMGetLabelByNum(dm, l, &label));
1225:     PetscCall(PetscObjectGetName((PetscObject)label, &name));
1226:     PetscCall(DMGetLabelOutput(dm, name, &output));
1227:     PetscCall(PetscStrncmp(name, "depth", 10, &isDepth));
1228:     if (isDepth || !output) continue;
1229:     PetscCall(PetscViewerHDF5PushGroup(viewer, name));
1230:     PetscCall(DMLabelGetValueIS(label, &valueIS));
1231:     /* Must copy to a new IS on the global comm */
1232:     PetscCall(ISGetLocalSize(valueIS, &numValues));
1233:     PetscCall(ISGetIndices(valueIS, &values));
1234:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS));
1235:     PetscCall(ISRestoreIndices(valueIS, &values));
1236:     PetscCall(ISAllGather(pvalueIS, &globalValueIS));
1237:     PetscCall(ISDestroy(&pvalueIS));
1238:     PetscCall(ISSortRemoveDups(globalValueIS));
1239:     PetscCall(ISGetLocalSize(globalValueIS, &numValues));
1240:     PetscCall(ISGetIndices(globalValueIS, &values));
1241:     for (v = 0; v < numValues; ++v) {
1242:       IS              stratumIS, globalStratumIS;
1243:       const PetscInt *spoints = NULL;
1244:       PetscInt       *gspoints, n = 0, gn, p;
1245:       const char     *iname = "indices";
1246:       char            group[PETSC_MAX_PATH_LEN];

1248:       PetscCall(PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v]));
1249:       PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1250:       PetscCall(DMLabelGetStratumIS(label, values[v], &stratumIS));

1252:       if (stratumIS) PetscCall(ISGetLocalSize(stratumIS, &n));
1253:       if (stratumIS) PetscCall(ISGetIndices(stratumIS, &spoints));
1254:       for (gn = 0, p = 0; p < n; ++p)
1255:         if (gpoint[spoints[p]] >= 0) ++gn;
1256:       PetscCall(PetscMalloc1(gn, &gspoints));
1257:       for (gn = 0, p = 0; p < n; ++p)
1258:         if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
1259:       if (stratumIS) PetscCall(ISRestoreIndices(stratumIS, &spoints));
1260:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS));
1261:       PetscCall(PetscObjectSetName((PetscObject)globalStratumIS, iname));

1263:       PetscCall(ISView(globalStratumIS, viewer));
1264:       PetscCall(ISDestroy(&globalStratumIS));
1265:       PetscCall(ISDestroy(&stratumIS));
1266:       PetscCall(PetscViewerHDF5PopGroup(viewer));
1267:     }
1268:     PetscCall(ISRestoreIndices(globalValueIS, &values));
1269:     PetscCall(ISDestroy(&globalValueIS));
1270:     PetscCall(ISDestroy(&valueIS));
1271:     PetscCall(PetscViewerHDF5PopGroup(viewer));
1272:   }
1273:   PetscCall(ISRestoreIndices(globalPointNumbers, &gpoint));
1274:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1275:   PetscFunctionReturn(PETSC_SUCCESS);
1276: }

1278: /* We only write cells and vertices. Does this screw up parallel reading? */
1279: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
1280: {
1281:   IS                globalPointNumbers;
1282:   PetscViewerFormat format;
1283:   PetscBool         viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE;

1285:   PetscFunctionBegin;
1286:   PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
1287:   PetscCall(DMPlexCoordinatesView_HDF5_Internal(dm, viewer));

1289:   PetscCall(PetscViewerGetFormat(viewer, &format));
1290:   switch (format) {
1291:   case PETSC_VIEWER_HDF5_VIZ:
1292:     viz_geom  = PETSC_TRUE;
1293:     xdmf_topo = PETSC_TRUE;
1294:     break;
1295:   case PETSC_VIEWER_HDF5_XDMF:
1296:     xdmf_topo = PETSC_TRUE;
1297:     break;
1298:   case PETSC_VIEWER_HDF5_PETSC:
1299:     petsc_topo = PETSC_TRUE;
1300:     break;
1301:   case PETSC_VIEWER_DEFAULT:
1302:   case PETSC_VIEWER_NATIVE:
1303:     viz_geom   = PETSC_TRUE;
1304:     xdmf_topo  = PETSC_TRUE;
1305:     petsc_topo = PETSC_TRUE;
1306:     break;
1307:   default:
1308:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
1309:   }

1311:   if (viz_geom) PetscCall(DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer));
1312:   if (xdmf_topo) PetscCall(DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer));
1313:   if (petsc_topo) {
1314:     PetscCall(DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer));
1315:     PetscCall(DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer));
1316:   }

1318:   PetscCall(ISDestroy(&globalPointNumbers));
1319:   PetscFunctionReturn(PETSC_SUCCESS);
1320: }

1322: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
1323: {
1324:   MPI_Comm     comm;
1325:   const char  *topologydm_name;
1326:   const char  *sectiondm_name;
1327:   PetscSection gsection;

1329:   PetscFunctionBegin;
1330:   PetscCall(PetscObjectGetComm((PetscObject)sectiondm, &comm));
1331:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1332:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
1333:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1334:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1335:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1336:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1337:   PetscCall(DMGetGlobalSection(sectiondm, &gsection));
1338:   /* Save raw section */
1339:   PetscCall(PetscSectionView(gsection, viewer));
1340:   /* Save plex wrapper */
1341:   {
1342:     PetscInt        pStart, pEnd, p, n;
1343:     IS              globalPointNumbers;
1344:     const PetscInt *gpoints;
1345:     IS              orderIS;
1346:     PetscInt       *order;

1348:     PetscCall(PetscSectionGetChart(gsection, &pStart, &pEnd));
1349:     PetscCall(DMPlexCreatePointNumbering(dm, &globalPointNumbers));
1350:     PetscCall(ISGetIndices(globalPointNumbers, &gpoints));
1351:     for (p = pStart, n = 0; p < pEnd; ++p)
1352:       if (gpoints[p] >= 0) n++;
1353:     /* "order" is an array of global point numbers.
1354:        When loading, it is used with topology/order array
1355:        to match section points with plex topology points. */
1356:     PetscCall(PetscMalloc1(n, &order));
1357:     for (p = pStart, n = 0; p < pEnd; ++p)
1358:       if (gpoints[p] >= 0) order[n++] = gpoints[p];
1359:     PetscCall(ISRestoreIndices(globalPointNumbers, &gpoints));
1360:     PetscCall(ISDestroy(&globalPointNumbers));
1361:     PetscCall(ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS));
1362:     PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
1363:     PetscCall(ISView(orderIS, viewer));
1364:     PetscCall(ISDestroy(&orderIS));
1365:   }
1366:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1367:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1368:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1369:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1370:   PetscFunctionReturn(PETSC_SUCCESS);
1371: }

1373: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1374: {
1375:   const char *topologydm_name;
1376:   const char *sectiondm_name;
1377:   const char *vec_name;
1378:   PetscInt    bs;

1380:   PetscFunctionBegin;
1381:   /* Check consistency */
1382:   {
1383:     PetscSF pointsf, pointsf1;

1385:     PetscCall(DMGetPointSF(dm, &pointsf));
1386:     PetscCall(DMGetPointSF(sectiondm, &pointsf1));
1387:     PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1388:   }
1389:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1390:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
1391:   PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
1392:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1393:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1394:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1395:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1396:   PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
1397:   PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
1398:   PetscCall(VecGetBlockSize(vec, &bs));
1399:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
1400:   PetscCall(VecSetBlockSize(vec, 1));
1401:   /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but,    */
1402:   /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view    */
1403:   /* is set to VecView_Plex, which would save vec in a predefined location.  */
1404:   /* To save vec in where we want, we create a new Vec (temp) with           */
1405:   /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1406:   {
1407:     Vec                temp;
1408:     const PetscScalar *array;
1409:     PetscLayout        map;

1411:     PetscCall(VecCreate(PetscObjectComm((PetscObject)vec), &temp));
1412:     PetscCall(PetscObjectSetName((PetscObject)temp, vec_name));
1413:     PetscCall(VecGetLayout(vec, &map));
1414:     PetscCall(VecSetLayout(temp, map));
1415:     PetscCall(VecSetUp(temp));
1416:     PetscCall(VecGetArrayRead(vec, &array));
1417:     PetscCall(VecPlaceArray(temp, array));
1418:     PetscCall(VecView(temp, viewer));
1419:     PetscCall(VecResetArray(temp));
1420:     PetscCall(VecRestoreArrayRead(vec, &array));
1421:     PetscCall(VecDestroy(&temp));
1422:   }
1423:   PetscCall(VecSetBlockSize(vec, bs));
1424:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1425:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1426:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1427:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1428:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1429:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1430:   PetscFunctionReturn(PETSC_SUCCESS);
1431: }

1433: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1434: {
1435:   MPI_Comm     comm;
1436:   const char  *topologydm_name;
1437:   const char  *sectiondm_name;
1438:   const char  *vec_name;
1439:   PetscSection section;
1440:   PetscBool    includesConstraints;
1441:   Vec          gvec;
1442:   PetscInt     m, bs;

1444:   PetscFunctionBegin;
1445:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1446:   /* Check consistency */
1447:   {
1448:     PetscSF pointsf, pointsf1;

1450:     PetscCall(DMGetPointSF(dm, &pointsf));
1451:     PetscCall(DMGetPointSF(sectiondm, &pointsf1));
1452:     PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1453:   }
1454:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1455:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
1456:   PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
1457:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
1458:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
1459:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
1460:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
1461:   PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
1462:   PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
1463:   PetscCall(VecGetBlockSize(vec, &bs));
1464:   PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs));
1465:   PetscCall(VecCreate(comm, &gvec));
1466:   PetscCall(PetscObjectSetName((PetscObject)gvec, vec_name));
1467:   PetscCall(DMGetGlobalSection(sectiondm, &section));
1468:   PetscCall(PetscSectionGetIncludesConstraints(section, &includesConstraints));
1469:   if (includesConstraints) PetscCall(PetscSectionGetStorageSize(section, &m));
1470:   else PetscCall(PetscSectionGetConstrainedStorageSize(section, &m));
1471:   PetscCall(VecSetSizes(gvec, m, PETSC_DECIDE));
1472:   PetscCall(VecSetUp(gvec));
1473:   PetscCall(DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec));
1474:   PetscCall(DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec));
1475:   PetscCall(VecView(gvec, viewer));
1476:   PetscCall(VecDestroy(&gvec));
1477:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1478:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1479:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1480:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1481:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1482:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1483:   PetscFunctionReturn(PETSC_SUCCESS);
1484: }

1486: struct _n_LoadLabelsCtx {
1487:   MPI_Comm    comm;
1488:   PetscMPIInt rank;
1489:   DM          dm;
1490:   PetscViewer viewer;
1491:   DMLabel     label;
1492:   PetscSF     sfXC;
1493:   PetscLayout layoutX;
1494: };
1495: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;

1497: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1498: {
1499:   PetscFunctionBegin;
1500:   PetscCall(PetscNew(ctx));
1501:   PetscCall(PetscObjectReference((PetscObject)((*ctx)->dm = dm)));
1502:   PetscCall(PetscObjectReference((PetscObject)((*ctx)->viewer = viewer)));
1503:   PetscCall(PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm));
1504:   PetscCallMPI(MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank));
1505:   (*ctx)->sfXC = sfXC;
1506:   if (sfXC) {
1507:     PetscInt nX;

1509:     PetscCall(PetscObjectReference((PetscObject)sfXC));
1510:     PetscCall(PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL));
1511:     PetscCall(PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX));
1512:   }
1513:   PetscFunctionReturn(PETSC_SUCCESS);
1514: }

1516: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1517: {
1518:   PetscFunctionBegin;
1519:   if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
1520:   PetscCall(DMDestroy(&(*ctx)->dm));
1521:   PetscCall(PetscViewerDestroy(&(*ctx)->viewer));
1522:   PetscCall(PetscSFDestroy(&(*ctx)->sfXC));
1523:   PetscCall(PetscLayoutDestroy(&(*ctx)->layoutX));
1524:   PetscCall(PetscFree(*ctx));
1525:   PetscFunctionReturn(PETSC_SUCCESS);
1526: }

1528: /*
1529:     A: on-disk points
1530:     X: global points [0, NX)
1531:     C: distributed plex points
1532: */
1533: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1534: {
1535:   MPI_Comm        comm    = ctx->comm;
1536:   PetscSF         sfXC    = ctx->sfXC;
1537:   PetscLayout     layoutX = ctx->layoutX;
1538:   PetscSF         sfXA;
1539:   const PetscInt *A_points;
1540:   PetscInt        nX, nC;
1541:   PetscInt        n;

1543:   PetscFunctionBegin;
1544:   PetscCall(PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL));
1545:   PetscCall(ISGetLocalSize(stratumIS, &n));
1546:   PetscCall(ISGetIndices(stratumIS, &A_points));
1547:   PetscCall(PetscSFCreate(comm, &sfXA));
1548:   PetscCall(PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points));
1549:   PetscCall(ISCreate(comm, newStratumIS));
1550:   PetscCall(ISSetType(*newStratumIS, ISGENERAL));
1551:   {
1552:     PetscInt   i;
1553:     PetscBool *A_mask, *X_mask, *C_mask;

1555:     PetscCall(PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask));
1556:     for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE;
1557:     PetscCall(PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE));
1558:     PetscCall(PetscSFReduceEnd(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE));
1559:     PetscCall(PetscSFBcastBegin(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR));
1560:     PetscCall(PetscSFBcastEnd(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR));
1561:     PetscCall(ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask));
1562:     PetscCall(PetscFree3(A_mask, X_mask, C_mask));
1563:   }
1564:   PetscCall(PetscSFDestroy(&sfXA));
1565:   PetscCall(ISRestoreIndices(stratumIS, &A_points));
1566:   PetscFunctionReturn(PETSC_SUCCESS);
1567: }

1569: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1570: {
1571:   LoadLabelsCtx   ctx    = (LoadLabelsCtx)op_data;
1572:   PetscViewer     viewer = ctx->viewer;
1573:   DMLabel         label  = ctx->label;
1574:   MPI_Comm        comm   = ctx->comm;
1575:   IS              stratumIS;
1576:   const PetscInt *ind;
1577:   PetscInt        value, N, i;

1579:   PetscCall(PetscOptionsStringToInt(vname, &value));
1580:   PetscCall(ISCreate(comm, &stratumIS));
1581:   PetscCall(PetscObjectSetName((PetscObject)stratumIS, "indices"));
1582:   PetscCall(PetscViewerHDF5PushGroup(viewer, vname)); /* labels/<lname>/<vname> */

1584:   if (!ctx->sfXC) {
1585:     /* Force serial load */
1586:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N));
1587:     PetscCall(PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0));
1588:     PetscCall(PetscLayoutSetSize(stratumIS->map, N));
1589:   }
1590:   PetscCall(ISLoad(stratumIS, viewer));

1592:   if (ctx->sfXC) {
1593:     IS newStratumIS;

1595:     PetscCallHDF5(ReadLabelStratumHDF5_Distribute_Private, (stratumIS, ctx, &newStratumIS));
1596:     PetscCall(ISDestroy(&stratumIS));
1597:     stratumIS = newStratumIS;
1598:   }

1600:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1601:   PetscCall(ISGetLocalSize(stratumIS, &N));
1602:   PetscCall(ISGetIndices(stratumIS, &ind));
1603:   for (i = 0; i < N; ++i) PetscCall(DMLabelSetValue(label, ind[i], value));
1604:   PetscCall(ISRestoreIndices(stratumIS, &ind));
1605:   PetscCall(ISDestroy(&stratumIS));
1606:   return 0;
1607: }

1609: /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */
1610: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1611: {
1612:   LoadLabelsCtx  ctx = (LoadLabelsCtx)op_data;
1613:   DM             dm  = ctx->dm;
1614:   hsize_t        idx = 0;
1615:   PetscErrorCode ierr;
1616:   PetscBool      flg;
1617:   herr_t         err;

1619:   PetscCall(DMHasLabel(dm, lname, &flg));
1620:   if (flg) PetscCall(DMRemoveLabel(dm, lname, NULL));
1621:   ierr = DMCreateLabel(dm, lname);
1622:   if (ierr) return (herr_t)ierr;
1623:   ierr = DMGetLabel(dm, lname, &ctx->label);
1624:   if (ierr) return (herr_t)ierr;
1625:   ierr = PetscViewerHDF5PushGroup(ctx->viewer, lname);
1626:   if (ierr) return (herr_t)ierr;
1627:   /* Iterate over the label's strata */
1628:   PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1629:   ierr = PetscViewerHDF5PopGroup(ctx->viewer);
1630:   if (ierr) return (herr_t)ierr;
1631:   return err;
1632: }

1634: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1635: {
1636:   const char          *topologydm_name;
1637:   LoadLabelsCtx        ctx;
1638:   hsize_t              idx = 0;
1639:   char                 group[PETSC_MAX_PATH_LEN];
1640:   DMPlexStorageVersion version;
1641:   PetscBool            distributed, hasGroup;

1643:   PetscFunctionBegin;
1644:   PetscCall(DMPlexIsDistributed(dm, &distributed));
1645:   if (distributed) PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
1646:   PetscCall(LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx));
1647:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
1648:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
1649:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1650:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name));
1651:   } else {
1652:     PetscCall(PetscStrncpy(group, "labels", sizeof(group)));
1653:   }
1654:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));
1655:   PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup));
1656:   if (hasGroup) {
1657:     hid_t fileId, groupId;

1659:     PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId));
1660:     /* Iterate over labels */
1661:     PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1662:     PetscCallHDF5(H5Gclose, (groupId));
1663:   }
1664:   PetscCall(PetscViewerHDF5PopGroup(viewer));
1665:   PetscCall(LoadLabelsCtxDestroy(&ctx));
1666:   PetscFunctionReturn(PETSC_SUCCESS);
1667: }

1669: static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm)
1670: {
1671:   MPI_Comm        comm;
1672:   PetscMPIInt     size, rank;
1673:   PetscInt        dist_size;
1674:   const char     *distribution_name;
1675:   PetscInt        p, lsize;
1676:   IS              chartSizesIS, ownersIS, gpointsIS;
1677:   const PetscInt *chartSize, *owners, *gpoints;
1678:   PetscLayout     layout;
1679:   PetscBool       has;

1681:   PetscFunctionBegin;
1682:   *distsf = NULL;
1683:   *distdm = NULL;
1684:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1685:   PetscCallMPI(MPI_Comm_size(comm, &size));
1686:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1687:   PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
1688:   if (!distribution_name) PetscFunctionReturn(PETSC_SUCCESS);
1689:   PetscCall(PetscLogEventBegin(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
1690:   PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &has));
1691:   if (!has) {
1692:     char *full_group;

1694:     PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &full_group));
1695:     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);
1696:   }
1697:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size));
1698:   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);
1699:   PetscCall(ISCreate(comm, &chartSizesIS));
1700:   PetscCall(PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes"));
1701:   PetscCall(ISCreate(comm, &ownersIS));
1702:   PetscCall(PetscObjectSetName((PetscObject)ownersIS, "owners"));
1703:   PetscCall(ISCreate(comm, &gpointsIS));
1704:   PetscCall(PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers"));
1705:   PetscCall(PetscLayoutSetLocalSize(chartSizesIS->map, 1));
1706:   PetscCall(ISLoad(chartSizesIS, viewer));
1707:   PetscCall(ISGetIndices(chartSizesIS, &chartSize));
1708:   PetscCall(PetscLayoutSetLocalSize(ownersIS->map, *chartSize));
1709:   PetscCall(PetscLayoutSetLocalSize(gpointsIS->map, *chartSize));
1710:   PetscCall(ISLoad(ownersIS, viewer));
1711:   PetscCall(ISLoad(gpointsIS, viewer));
1712:   PetscCall(ISGetIndices(ownersIS, &owners));
1713:   PetscCall(ISGetIndices(gpointsIS, &gpoints));
1714:   PetscCall(PetscSFCreate(comm, distsf));
1715:   PetscCall(PetscSFSetFromOptions(*distsf));
1716:   PetscCall(PetscLayoutCreate(comm, &layout));
1717:   PetscCall(PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL));
1718:   PetscCall(PetscLayoutSetLocalSize(layout, lsize));
1719:   PetscCall(PetscLayoutSetBlockSize(layout, 1));
1720:   PetscCall(PetscLayoutSetUp(layout));
1721:   PetscCall(PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints));
1722:   PetscCall(PetscLayoutDestroy(&layout));
1723:   /* Migrate DM */
1724:   {
1725:     PetscInt     pStart, pEnd;
1726:     PetscSFNode *buffer0, *buffer1, *buffer2;

1728:     PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1729:     PetscCall(PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1));
1730:     PetscCall(PetscMalloc1(*chartSize, &buffer2));
1731:     {
1732:       PetscSF            workPointSF;
1733:       PetscInt           workNroots, workNleaves;
1734:       const PetscInt    *workIlocal;
1735:       const PetscSFNode *workIremote;

1737:       for (p = pStart; p < pEnd; ++p) {
1738:         buffer0[p - pStart].rank  = rank;
1739:         buffer0[p - pStart].index = p - pStart;
1740:       }
1741:       PetscCall(DMGetPointSF(dm, &workPointSF));
1742:       PetscCall(PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote));
1743:       for (p = 0; p < workNleaves; ++p) {
1744:         PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p);

1746:         buffer0[workIlocalp].rank = -1;
1747:       }
1748:     }
1749:     for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1750:     for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1;
1751:     PetscCall(PetscSFReduceBegin(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1752:     PetscCall(PetscSFReduceEnd(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1753:     PetscCall(PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE));
1754:     PetscCall(PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE));
1755:     if (PetscDefined(USE_DEBUG)) {
1756:       for (p = 0; p < *chartSize; ++p) {
1757:         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);
1758:       }
1759:     }
1760:     PetscCall(PetscFree2(buffer0, buffer1));
1761:     PetscCall(DMCreate(comm, distdm));
1762:     PetscCall(DMSetType(*distdm, DMPLEX));
1763:     PetscCall(PetscObjectSetName((PetscObject)(*distdm), ((PetscObject)dm)->name));
1764:     PetscCall(DMPlexDistributionSetName(*distdm, distribution_name));
1765:     {
1766:       PetscSF migrationSF;

1768:       PetscCall(PetscSFCreate(comm, &migrationSF));
1769:       PetscCall(PetscSFSetFromOptions(migrationSF));
1770:       PetscCall(PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER));
1771:       PetscCall(PetscSFSetUp(migrationSF));
1772:       PetscCall(DMPlexMigrate(dm, migrationSF, *distdm));
1773:       PetscCall(PetscSFDestroy(&migrationSF));
1774:     }
1775:   }
1776:   /* Set pointSF */
1777:   {
1778:     PetscSF      pointSF;
1779:     PetscInt    *ilocal, nleaves, q;
1780:     PetscSFNode *iremote, *buffer0, *buffer1;

1782:     PetscCall(PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1));
1783:     for (p = 0, nleaves = 0; p < *chartSize; ++p) {
1784:       if (owners[p] == rank) {
1785:         buffer0[p].rank = rank;
1786:       } else {
1787:         buffer0[p].rank = -1;
1788:         nleaves++;
1789:       }
1790:       buffer0[p].index = p;
1791:     }
1792:     for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1793:     PetscCall(PetscSFReduceBegin(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1794:     PetscCall(PetscSFReduceEnd(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC));
1795:     for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1;
1796:     PetscCall(PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE));
1797:     PetscCall(PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE));
1798:     if (PetscDefined(USE_DEBUG)) {
1799:       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);
1800:     }
1801:     PetscCall(PetscMalloc1(nleaves, &ilocal));
1802:     PetscCall(PetscMalloc1(nleaves, &iremote));
1803:     for (p = 0, q = 0; p < *chartSize; ++p) {
1804:       if (buffer0[p].rank != rank) {
1805:         ilocal[q]        = p;
1806:         iremote[q].rank  = buffer0[p].rank;
1807:         iremote[q].index = buffer0[p].index;
1808:         q++;
1809:       }
1810:     }
1811:     PetscCheck(q == nleaves, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching leaf sizes: %" PetscInt_FMT " != %" PetscInt_FMT, q, nleaves);
1812:     PetscCall(PetscFree2(buffer0, buffer1));
1813:     PetscCall(PetscSFCreate(comm, &pointSF));
1814:     PetscCall(PetscSFSetFromOptions(pointSF));
1815:     PetscCall(PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
1816:     PetscCall(DMSetPointSF(*distdm, pointSF));
1817:     {
1818:       DM cdm;

1820:       PetscCall(DMGetCoordinateDM(*distdm, &cdm));
1821:       PetscCall(DMSetPointSF(cdm, pointSF));
1822:     }
1823:     PetscCall(PetscSFDestroy(&pointSF));
1824:   }
1825:   PetscCall(ISRestoreIndices(chartSizesIS, &chartSize));
1826:   PetscCall(ISRestoreIndices(ownersIS, &owners));
1827:   PetscCall(ISRestoreIndices(gpointsIS, &gpoints));
1828:   PetscCall(ISDestroy(&chartSizesIS));
1829:   PetscCall(ISDestroy(&ownersIS));
1830:   PetscCall(ISDestroy(&gpointsIS));
1831:   /* Record that overlap has been manually created.               */
1832:   /* This is to pass `DMPlexCheckPointSF()`, which checks that    */
1833:   /* pointSF does not contain cells in the leaves if overlap = 0. */
1834:   PetscCall(DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL));
1835:   PetscCall(DMPlexDistributeSetDefault(*distdm, PETSC_FALSE));
1836:   PetscCall(DMPlexReorderSetDefault(*distdm, DMPLEX_REORDER_DEFAULT_FALSE));
1837:   PetscCall(PetscLogEventEnd(DMPLEX_DistributionLoad, viewer, 0, 0, 0));
1838:   PetscFunctionReturn(PETSC_SUCCESS);
1839: }

1841: static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf)
1842: {
1843:   MPI_Comm        comm;
1844:   const char     *pointsName, *coneSizesName, *conesName, *orientationsName;
1845:   IS              pointsIS, coneSizesIS, conesIS, orientationsIS;
1846:   const PetscInt *points, *coneSizes, *cones, *orientations;
1847:   PetscInt       *cone, *ornt;
1848:   PetscInt        dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1849:   PetscMPIInt     size, rank;

1851:   PetscFunctionBegin;
1852:   pointsName       = "order";
1853:   coneSizesName    = "cones";
1854:   conesName        = "cells";
1855:   orientationsName = "orientation";
1856:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1857:   PetscCallMPI(MPI_Comm_size(comm, &size));
1858:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1859:   PetscCall(ISCreate(comm, &pointsIS));
1860:   PetscCall(PetscObjectSetName((PetscObject)pointsIS, pointsName));
1861:   PetscCall(ISCreate(comm, &coneSizesIS));
1862:   PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
1863:   PetscCall(ISCreate(comm, &conesIS));
1864:   PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
1865:   PetscCall(ISCreate(comm, &orientationsIS));
1866:   PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
1867:   PetscCall(PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim));
1868:   PetscCall(DMSetDimension(dm, dim));
1869:   {
1870:     /* Force serial load */
1871:     PetscCall(PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np));
1872:     PetscCall(PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0));
1873:     PetscCall(PetscLayoutSetSize(pointsIS->map, Np));
1874:     pEnd = rank == 0 ? Np : 0;
1875:     PetscCall(PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np));
1876:     PetscCall(PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0));
1877:     PetscCall(PetscLayoutSetSize(coneSizesIS->map, Np));
1878:     PetscCall(PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N));
1879:     PetscCall(PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0));
1880:     PetscCall(PetscLayoutSetSize(conesIS->map, N));
1881:     PetscCall(PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N));
1882:     PetscCall(PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0));
1883:     PetscCall(PetscLayoutSetSize(orientationsIS->map, N));
1884:   }
1885:   PetscCall(ISLoad(pointsIS, viewer));
1886:   PetscCall(ISLoad(coneSizesIS, viewer));
1887:   PetscCall(ISLoad(conesIS, viewer));
1888:   PetscCall(ISLoad(orientationsIS, viewer));
1889:   /* Create Plex */
1890:   PetscCall(DMPlexSetChart(dm, 0, pEnd));
1891:   PetscCall(ISGetIndices(pointsIS, &points));
1892:   PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
1893:   for (p = 0; p < pEnd; ++p) {
1894:     PetscCall(DMPlexSetConeSize(dm, points[p], coneSizes[p]));
1895:     maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
1896:   }
1897:   PetscCall(DMSetUp(dm));
1898:   PetscCall(ISGetIndices(conesIS, &cones));
1899:   PetscCall(ISGetIndices(orientationsIS, &orientations));
1900:   PetscCall(PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt));
1901:   for (p = 0, q = 0; p < pEnd; ++p) {
1902:     for (c = 0; c < coneSizes[p]; ++c, ++q) {
1903:       cone[c] = cones[q];
1904:       ornt[c] = orientations[q];
1905:     }
1906:     PetscCall(DMPlexSetCone(dm, points[p], cone));
1907:     PetscCall(DMPlexSetConeOrientation(dm, points[p], ornt));
1908:   }
1909:   PetscCall(PetscFree2(cone, ornt));
1910:   /* Create global section migration SF */
1911:   if (sf) {
1912:     PetscLayout layout;
1913:     PetscInt   *globalIndices;

1915:     PetscCall(PetscMalloc1(pEnd, &globalIndices));
1916:     /* plex point == globalPointNumber in this case */
1917:     for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1918:     PetscCall(PetscLayoutCreate(comm, &layout));
1919:     PetscCall(PetscLayoutSetSize(layout, Np));
1920:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
1921:     PetscCall(PetscLayoutSetUp(layout));
1922:     PetscCall(PetscSFCreate(comm, sf));
1923:     PetscCall(PetscSFSetFromOptions(*sf));
1924:     PetscCall(PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices));
1925:     PetscCall(PetscLayoutDestroy(&layout));
1926:     PetscCall(PetscFree(globalIndices));
1927:   }
1928:   /* Clean-up */
1929:   PetscCall(ISRestoreIndices(pointsIS, &points));
1930:   PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
1931:   PetscCall(ISRestoreIndices(conesIS, &cones));
1932:   PetscCall(ISRestoreIndices(orientationsIS, &orientations));
1933:   PetscCall(ISDestroy(&pointsIS));
1934:   PetscCall(ISDestroy(&coneSizesIS));
1935:   PetscCall(ISDestroy(&conesIS));
1936:   PetscCall(ISDestroy(&orientationsIS));
1937:   /* Fill in the rest of the topology structure */
1938:   PetscCall(DMPlexSymmetrize(dm));
1939:   PetscCall(DMPlexStratify(dm));
1940:   PetscFunctionReturn(PETSC_SUCCESS);
1941: }

1943: /* Representation of two DMPlex strata in 0-based global numbering */
1944: struct _n_PlexLayer {
1945:   PetscInt     d;
1946:   IS           conesIS, orientationsIS;
1947:   PetscSection coneSizesSection;
1948:   PetscLayout  vertexLayout;
1949:   PetscSF      overlapSF, l2gSF; //TODO maybe confusing names (in DMPlex in general)
1950:   PetscInt     offset, conesOffset, leafOffset;
1951: };
1952: typedef struct _n_PlexLayer *PlexLayer;

1954: static PetscErrorCode PlexLayerDestroy(PlexLayer *layer)
1955: {
1956:   PetscFunctionBegin;
1957:   if (!*layer) PetscFunctionReturn(PETSC_SUCCESS);
1958:   PetscCall(PetscSectionDestroy(&(*layer)->coneSizesSection));
1959:   PetscCall(ISDestroy(&(*layer)->conesIS));
1960:   PetscCall(ISDestroy(&(*layer)->orientationsIS));
1961:   PetscCall(PetscSFDestroy(&(*layer)->overlapSF));
1962:   PetscCall(PetscSFDestroy(&(*layer)->l2gSF));
1963:   PetscCall(PetscLayoutDestroy(&(*layer)->vertexLayout));
1964:   PetscCall(PetscFree(*layer));
1965:   PetscFunctionReturn(PETSC_SUCCESS);
1966: }

1968: static PetscErrorCode PlexLayerCreate_Private(PlexLayer *layer)
1969: {
1970:   PetscFunctionBegin;
1971:   PetscCall(PetscNew(layer));
1972:   (*layer)->d           = -1;
1973:   (*layer)->offset      = -1;
1974:   (*layer)->conesOffset = -1;
1975:   (*layer)->leafOffset  = -1;
1976:   PetscFunctionReturn(PETSC_SUCCESS);
1977: }

1979: static PetscErrorCode PlexLayerLoad_Private(PlexLayer layer, PetscViewer viewer, PetscInt d, PetscLayout pointsLayout)
1980: {
1981:   char         path[128];
1982:   MPI_Comm     comm;
1983:   const char  *coneSizesName, *conesName, *orientationsName;
1984:   IS           coneSizesIS, conesIS, orientationsIS;
1985:   PetscSection coneSizesSection;
1986:   PetscLayout  vertexLayout = NULL;
1987:   PetscInt     s;

1989:   PetscFunctionBegin;
1990:   coneSizesName    = "cone_sizes";
1991:   conesName        = "cones";
1992:   orientationsName = "orientations";
1993:   PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));

1995:   /* query size of next lower depth stratum (next lower dimension) */
1996:   if (d > 0) {
1997:     PetscInt NVertices;

1999:     PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT "/%s", d - 1, coneSizesName));
2000:     PetscCall(PetscViewerHDF5ReadSizes(viewer, path, NULL, &NVertices));
2001:     PetscCall(PetscLayoutCreate(comm, &vertexLayout));
2002:     PetscCall(PetscLayoutSetSize(vertexLayout, NVertices));
2003:     PetscCall(PetscLayoutSetUp(vertexLayout));
2004:   }

2006:   PetscCall(PetscSNPrintf(path, sizeof(path), "%" PetscInt_FMT, d));
2007:   PetscCall(PetscViewerHDF5PushGroup(viewer, path));

2009:   /* create coneSizesSection from stored IS coneSizes */
2010:   {
2011:     const PetscInt *coneSizes;

2013:     PetscCall(ISCreate(comm, &coneSizesIS));
2014:     PetscCall(PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName));
2015:     if (pointsLayout) PetscCall(ISSetLayout(coneSizesIS, pointsLayout));
2016:     PetscCall(ISLoad(coneSizesIS, viewer));
2017:     if (!pointsLayout) PetscCall(ISGetLayout(coneSizesIS, &pointsLayout));
2018:     PetscCall(ISGetIndices(coneSizesIS, &coneSizes));
2019:     PetscCall(PetscSectionCreate(comm, &coneSizesSection));
2020:     //TODO different start ?
2021:     PetscCall(PetscSectionSetChart(coneSizesSection, 0, pointsLayout->n));
2022:     for (s = 0; s < pointsLayout->n; ++s) PetscCall(PetscSectionSetDof(coneSizesSection, s, coneSizes[s]));
2023:     PetscCall(PetscSectionSetUp(coneSizesSection));
2024:     PetscCall(ISRestoreIndices(coneSizesIS, &coneSizes));
2025:     {
2026:       PetscLayout tmp = NULL;

2028:       /* We need to keep the layout until the end of function */
2029:       PetscCall(PetscLayoutReference((PetscLayout)pointsLayout, &tmp));
2030:     }
2031:     PetscCall(ISDestroy(&coneSizesIS));
2032:   }

2034:   /* use value layout of coneSizesSection as layout of cones and orientations */
2035:   {
2036:     PetscLayout conesLayout;

2038:     PetscCall(PetscSectionGetValueLayout(comm, coneSizesSection, &conesLayout));
2039:     PetscCall(ISCreate(comm, &conesIS));
2040:     PetscCall(ISCreate(comm, &orientationsIS));
2041:     PetscCall(PetscObjectSetName((PetscObject)conesIS, conesName));
2042:     PetscCall(PetscObjectSetName((PetscObject)orientationsIS, orientationsName));
2043:     PetscCall(PetscLayoutDuplicate(conesLayout, &conesIS->map));
2044:     PetscCall(PetscLayoutDuplicate(conesLayout, &orientationsIS->map));
2045:     PetscCall(ISLoad(conesIS, viewer));
2046:     PetscCall(ISLoad(orientationsIS, viewer));
2047:     PetscCall(PetscLayoutDestroy(&conesLayout));
2048:   }

2050:   /* check assertion that layout of points is the same as point layout of coneSizesSection */
2051:   {
2052:     PetscLayout pointsLayout0;
2053:     PetscBool   flg;

2055:     PetscCall(PetscSectionGetPointLayout(comm, coneSizesSection, &pointsLayout0));
2056:     PetscCall(PetscLayoutCompare(pointsLayout, pointsLayout0, &flg));
2057:     PetscCheck(flg, comm, PETSC_ERR_PLIB, "points layout != coneSizesSection point layout");
2058:     PetscCall(PetscLayoutDestroy(&pointsLayout0));
2059:   }
2060:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2061:   PetscCall(PetscLayoutDestroy(&pointsLayout));

2063:   layer->d                = d;
2064:   layer->conesIS          = conesIS;
2065:   layer->coneSizesSection = coneSizesSection;
2066:   layer->orientationsIS   = orientationsIS;
2067:   layer->vertexLayout     = vertexLayout;
2068:   PetscFunctionReturn(PETSC_SUCCESS);
2069: }

2071: static PetscErrorCode PlexLayerDistribute_Private(PlexLayer layer, PetscSF cellLocalToGlobalSF)
2072: {
2073:   IS           newConesIS, newOrientationsIS;
2074:   PetscSection newConeSizesSection;
2075:   MPI_Comm     comm;

2077:   PetscFunctionBegin;
2078:   PetscCall(PetscObjectGetComm((PetscObject)cellLocalToGlobalSF, &comm));
2079:   PetscCall(PetscSectionCreate(comm, &newConeSizesSection));
2080:   //TODO rename to something like ISDistribute() with optional PetscSection argument
2081:   PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->conesIS, newConeSizesSection, &newConesIS));
2082:   PetscCall(DMPlexDistributeFieldIS(NULL, cellLocalToGlobalSF, layer->coneSizesSection, layer->orientationsIS, newConeSizesSection, &newOrientationsIS));

2084:   PetscCall(PetscObjectSetName((PetscObject)newConeSizesSection, ((PetscObject)layer->coneSizesSection)->name));
2085:   PetscCall(PetscObjectSetName((PetscObject)newConesIS, ((PetscObject)layer->conesIS)->name));
2086:   PetscCall(PetscObjectSetName((PetscObject)newOrientationsIS, ((PetscObject)layer->orientationsIS)->name));
2087:   PetscCall(PetscSectionDestroy(&layer->coneSizesSection));
2088:   PetscCall(ISDestroy(&layer->conesIS));
2089:   PetscCall(ISDestroy(&layer->orientationsIS));
2090:   layer->coneSizesSection = newConeSizesSection;
2091:   layer->conesIS          = newConesIS;
2092:   layer->orientationsIS   = newOrientationsIS;
2093:   PetscFunctionReturn(PETSC_SUCCESS);
2094: }

2096:   //TODO share code with DMPlexBuildFromCellListParallel()
2097: #include <petsc/private/hashseti.h>
2098: static PetscErrorCode PlexLayerCreateSFs_Private(PlexLayer layer, PetscSF *vertexOverlapSF, PetscSF *sfXC)
2099: {
2100:   PetscLayout     vertexLayout     = layer->vertexLayout;
2101:   PetscSection    coneSection      = layer->coneSizesSection;
2102:   IS              cellVertexData   = layer->conesIS;
2103:   IS              coneOrientations = layer->orientationsIS;
2104:   PetscSF         vl2gSF, vOverlapSF;
2105:   PetscInt       *verticesAdj;
2106:   PetscInt        i, n, numVerticesAdj;
2107:   const PetscInt *cvd, *co = NULL;
2108:   MPI_Comm        comm;

2110:   PetscFunctionBegin;
2111:   PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
2112:   PetscCall(PetscSectionGetStorageSize(coneSection, &n));
2113:   {
2114:     PetscInt n0;

2116:     PetscCall(ISGetLocalSize(cellVertexData, &n0));
2117:     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);
2118:     PetscCall(ISGetIndices(cellVertexData, &cvd));
2119:   }
2120:   if (coneOrientations) {
2121:     PetscInt n0;

2123:     PetscCall(ISGetLocalSize(coneOrientations, &n0));
2124:     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);
2125:     PetscCall(ISGetIndices(coneOrientations, &co));
2126:   }
2127:   /* Get/check global number of vertices */
2128:   {
2129:     PetscInt NVerticesInCells = PETSC_MIN_INT;

2131:     /* NVerticesInCells = max(cellVertexData) + 1 */
2132:     for (i = 0; i < n; i++)
2133:       if (cvd[i] > NVerticesInCells) NVerticesInCells = cvd[i];
2134:     ++NVerticesInCells;
2135:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, comm));

2137:     if (vertexLayout->n == PETSC_DECIDE && vertexLayout->N == PETSC_DECIDE) vertexLayout->N = NVerticesInCells;
2138:     else
2139:       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,
2140:                  vertexLayout->N, NVerticesInCells);
2141:     PetscCall(PetscLayoutSetUp(vertexLayout));
2142:   }
2143:   /* Find locally unique vertices in cellVertexData */
2144:   {
2145:     PetscHSetI vhash;
2146:     PetscInt   off = 0;

2148:     PetscCall(PetscHSetICreate(&vhash));
2149:     for (i = 0; i < n; ++i) PetscCall(PetscHSetIAdd(vhash, cvd[i]));
2150:     PetscCall(PetscHSetIGetSize(vhash, &numVerticesAdj));
2151:     PetscCall(PetscMalloc1(numVerticesAdj, &verticesAdj));
2152:     PetscCall(PetscHSetIGetElems(vhash, &off, verticesAdj));
2153:     PetscAssert(off == numVerticesAdj, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assertion failed: off == numVerticesAdj (%" PetscInt_FMT " != %" PetscInt_FMT ")", off, numVerticesAdj);
2154:     PetscCall(PetscHSetIDestroy(&vhash));
2155:   }
2156:   /* We must sort vertices to preserve numbering */
2157:   PetscCall(PetscSortInt(numVerticesAdj, verticesAdj));
2158:   /* Connect locally unique vertices with star forest */
2159:   PetscCall(PetscSFCreateByMatchingIndices(vertexLayout, numVerticesAdj, verticesAdj, NULL, 0, numVerticesAdj, verticesAdj, NULL, 0, &vl2gSF, &vOverlapSF));
2160:   PetscCall(PetscObjectSetName((PetscObject)vOverlapSF, "overlapSF"));
2161:   PetscCall(PetscObjectSetName((PetscObject)vl2gSF, "localToGlobalSF"));

2163:   PetscCall(PetscFree(verticesAdj));
2164:   *vertexOverlapSF = vOverlapSF;
2165:   *sfXC            = vl2gSF;
2166:   PetscFunctionReturn(PETSC_SUCCESS);
2167: }

2169: static PetscErrorCode PlexLayerCreateCellSFs_Private(PlexLayer layer, PetscSF *cellOverlapSF, PetscSF *cellLocalToGlobalSF)
2170: {
2171:   PetscSection coneSection = layer->coneSizesSection;
2172:   PetscInt     nCells;
2173:   MPI_Comm     comm;

2175:   PetscFunctionBegin;
2176:   PetscCall(PetscObjectGetComm((PetscObject)coneSection, &comm));
2177:   {
2178:     PetscInt cStart;

2180:     PetscCall(PetscSectionGetChart(coneSection, &cStart, &nCells));
2181:     PetscCheck(cStart == 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "coneSection must start at 0");
2182:   }
2183:   /* Create overlapSF as empty SF with the right number of roots */
2184:   PetscCall(PetscSFCreate(comm, cellOverlapSF));
2185:   PetscCall(PetscSFSetGraph(*cellOverlapSF, nCells, 0, NULL, PETSC_USE_POINTER, NULL, PETSC_USE_POINTER));
2186:   PetscCall(PetscSFSetUp(*cellOverlapSF));
2187:   /* Create localToGlobalSF as identity mapping */
2188:   {
2189:     PetscLayout map;

2191:     PetscCall(PetscLayoutCreateFromSizes(comm, nCells, PETSC_DECIDE, 1, &map));
2192:     PetscCall(PetscSFCreateFromLayouts(map, map, cellLocalToGlobalSF));
2193:     PetscCall(PetscSFSetUp(*cellLocalToGlobalSF));
2194:     PetscCall(PetscLayoutDestroy(&map));
2195:   }
2196:   PetscFunctionReturn(PETSC_SUCCESS);
2197: }

2199: static PetscErrorCode DMPlexTopologyBuildFromLayers_Private(DM dm, PetscInt depth, PlexLayer *layers, IS strataPermutation)
2200: {
2201:   const PetscInt *permArr;
2202:   PetscInt        d, nPoints;
2203:   MPI_Comm        comm;

2205:   PetscFunctionBegin;
2206:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2207:   PetscCall(ISGetIndices(strataPermutation, &permArr));

2209:   /* Count points, strata offsets and cones offsets (taking strataPermutation into account) */
2210:   {
2211:     PetscInt stratumOffset = 0;
2212:     PetscInt conesOffset   = 0;

2214:     for (d = 0; d <= depth; d++) {
2215:       const PetscInt  e = permArr[d];
2216:       const PlexLayer l = layers[e];
2217:       PetscInt        lo, n, size;

2219:       PetscCall(PetscSectionGetChart(l->coneSizesSection, &lo, &n));
2220:       PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &size));
2221:       PetscCheck(lo == 0, comm, PETSC_ERR_PLIB, "starting point should be 0 in coneSizesSection %" PetscInt_FMT, d);
2222:       l->offset      = stratumOffset;
2223:       l->conesOffset = conesOffset;
2224:       stratumOffset += n;
2225:       conesOffset += size;
2226:     }
2227:     nPoints = stratumOffset;
2228:   }

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

2234:   /* Set up plex coneSection from layer coneSections */
2235:   {
2236:     PetscSection coneSection;

2238:     PetscCall(DMPlexGetConeSection(dm, &coneSection));
2239:     for (d = 0; d <= depth; d++) {
2240:       const PlexLayer l = layers[d];
2241:       PetscInt        n, q;

2243:       PetscCall(PetscSectionGetChart(l->coneSizesSection, NULL, &n));
2244:       for (q = 0; q < n; q++) {
2245:         const PetscInt p = l->offset + q;
2246:         PetscInt       coneSize;

2248:         PetscCall(PetscSectionGetDof(l->coneSizesSection, q, &coneSize));
2249:         PetscCall(PetscSectionSetDof(coneSection, p, coneSize));
2250:       }
2251:     }
2252:   }
2253:   //TODO this is terrible, DMSetUp_Plex() should be DMPlexSetUpSections() or so
2254:   PetscCall(DMSetUp(dm));

2256:   /* Renumber cones points from layer-global numbering to plex-local numbering */
2257:   {
2258:     PetscInt *cones, *ornts;

2260:     PetscCall(DMPlexGetCones(dm, &cones));
2261:     PetscCall(DMPlexGetConeOrientations(dm, &ornts));
2262:     for (d = 1; d <= depth; d++) {
2263:       const PlexLayer l = layers[d];
2264:       PetscInt        i, lConesSize;
2265:       PetscInt       *lCones;
2266:       const PetscInt *lOrnts;
2267:       PetscInt       *pCones = &cones[l->conesOffset];
2268:       PetscInt       *pOrnts = &ornts[l->conesOffset];

2270:       PetscCall(PetscSectionGetStorageSize(l->coneSizesSection, &lConesSize));
2271:       /* Get cones in local plex numbering */
2272:       {
2273:         ISLocalToGlobalMapping l2g;
2274:         PetscLayout            vertexLayout = l->vertexLayout;
2275:         PetscSF                vertexSF     = layers[d - 1]->l2gSF; /* vertices of this layer are cells of previous layer */
2276:         const PetscInt        *gCones;
2277:         PetscInt               lConesSize0;

2279:         PetscCall(ISGetLocalSize(l->conesIS, &lConesSize0));
2280:         PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(conesIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);
2281:         PetscCall(ISGetLocalSize(l->orientationsIS, &lConesSize0));
2282:         PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "layer %" PetscInt_FMT " size(orientationsIS) = %" PetscInt_FMT " != %" PetscInt_FMT " = storageSize(coneSizesSection)", d, lConesSize0, lConesSize);

2284:         PetscCall(PetscMalloc1(lConesSize, &lCones));
2285:         PetscCall(ISGetIndices(l->conesIS, &gCones));
2286:         PetscCall(ISLocalToGlobalMappingCreateSF(vertexSF, vertexLayout->rstart, &l2g));
2287:         PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, lConesSize, gCones, &lConesSize0, lCones));
2288:         PetscCheck(lConesSize0 == lConesSize, comm, PETSC_ERR_PLIB, "global to local does not cover all indices (%" PetscInt_FMT " of %" PetscInt_FMT ")", lConesSize0, lConesSize);
2289:         PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
2290:         PetscCall(ISRestoreIndices(l->conesIS, &gCones));
2291:       }
2292:       PetscCall(ISGetIndices(l->orientationsIS, &lOrnts));
2293:       /* Set cones, need to add stratum offset */
2294:       for (i = 0; i < lConesSize; i++) {
2295:         pCones[i] = lCones[i] + layers[d - 1]->offset; /* cone points of current layer are points of previous layer */
2296:         pOrnts[i] = lOrnts[i];
2297:       }
2298:       PetscCall(PetscFree(lCones));
2299:       PetscCall(ISRestoreIndices(l->orientationsIS, &lOrnts));
2300:     }
2301:   }
2302:   PetscCall(DMPlexSymmetrize(dm));
2303:   PetscCall(DMPlexStratify(dm));
2304:   PetscCall(ISRestoreIndices(strataPermutation, &permArr));
2305:   PetscFunctionReturn(PETSC_SUCCESS);
2306: }

2308: static PetscErrorCode PlexLayerConcatenateSFs_Private(MPI_Comm comm, PetscInt depth, PlexLayer layers[], IS strataPermutation, PetscSF *overlapSF, PetscSF *l2gSF)
2309: {
2310:   PetscInt        d;
2311:   PetscSF        *osfs, *lsfs;
2312:   PetscInt       *leafOffsets;
2313:   const PetscInt *permArr;

2315:   PetscFunctionBegin;
2316:   PetscCall(ISGetIndices(strataPermutation, &permArr));
2317:   PetscCall(PetscCalloc3(depth + 1, &osfs, depth + 1, &lsfs, depth + 1, &leafOffsets));
2318:   for (d = 0; d <= depth; d++) {
2319:     const PetscInt e = permArr[d];

2321:     PetscAssert(e == layers[e]->d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "assertion: e == layers[e]->d");
2322:     osfs[d]        = layers[e]->overlapSF;
2323:     lsfs[d]        = layers[e]->l2gSF;
2324:     leafOffsets[d] = layers[e]->offset;
2325:   }
2326:   PetscCall(PetscSFConcatenate(comm, depth + 1, osfs, PETSCSF_CONCATENATE_../../../..MODE_LOCAL, leafOffsets, overlapSF));
2327:   PetscCall(PetscSFConcatenate(comm, depth + 1, lsfs, PETSCSF_CONCATENATE_../../../..MODE_GLOBAL, leafOffsets, l2gSF));
2328:   PetscCall(PetscFree3(osfs, lsfs, leafOffsets));
2329:   PetscCall(ISRestoreIndices(strataPermutation, &permArr));
2330:   PetscFunctionReturn(PETSC_SUCCESS);
2331: }

2333: static PetscErrorCode DMPlexTopologyLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF *sfXC)
2334: {
2335:   PlexLayer  *layers;
2336:   IS          strataPermutation;
2337:   PetscLayout pointsLayout;
2338:   PetscInt    depth;
2339:   PetscInt    d;
2340:   MPI_Comm    comm;

2342:   PetscFunctionBegin;
2343:   {
2344:     PetscInt dim;

2346:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "depth", PETSC_INT, NULL, &depth));
2347:     PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "cell_dim", PETSC_INT, NULL, &dim));
2348:     PetscCall(DMSetDimension(dm, dim));
2349:   }
2350:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));

2352:   {
2353:     IS spOnComm;

2355:     PetscCall(ISCreate(comm, &spOnComm));
2356:     PetscCall(PetscObjectSetName((PetscObject)spOnComm, "permutation"));
2357:     PetscCall(ISLoad(spOnComm, viewer));
2358:     /* have the same serial IS on every rank */
2359:     PetscCall(ISAllGather(spOnComm, &strataPermutation));
2360:     PetscCall(PetscObjectSetName((PetscObject)strataPermutation, ((PetscObject)spOnComm)->name));
2361:     PetscCall(ISDestroy(&spOnComm));
2362:   }

2364:   /* Create layers, load raw data for each layer */
2365:   PetscCall(PetscViewerHDF5PushGroup(viewer, "strata"));
2366:   PetscCall(PetscMalloc1(depth + 1, &layers));
2367:   for (d = depth, pointsLayout = NULL; d >= 0; pointsLayout = layers[d]->vertexLayout, d--) {
2368:     PetscCall(PlexLayerCreate_Private(&layers[d]));
2369:     PetscCall(PlexLayerLoad_Private(layers[d], viewer, d, pointsLayout));
2370:   }
2371:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* strata */

2373:   for (d = depth; d >= 0; d--) {
2374:     /* Redistribute cells and vertices for each applicable layer */
2375:     if (d < depth) PetscCall(PlexLayerDistribute_Private(layers[d], layers[d]->l2gSF));
2376:     /* Create vertex overlap SF and vertex localToGlobal SF for each applicable layer */
2377:     if (d > 0) PetscCall(PlexLayerCreateSFs_Private(layers[d], &layers[d - 1]->overlapSF, &layers[d - 1]->l2gSF));
2378:   }
2379:   /* Build trivial SFs for the cell layer as well */
2380:   PetscCall(PlexLayerCreateCellSFs_Private(layers[depth], &layers[depth]->overlapSF, &layers[depth]->l2gSF));

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

2385:   /* Build overall point SF alias overlap SF */
2386:   {
2387:     PetscSF overlapSF;

2389:     PetscCall(PlexLayerConcatenateSFs_Private(comm, depth, layers, strataPermutation, &overlapSF, sfXC));
2390:     PetscCall(DMSetPointSF(dm, overlapSF));
2391:     PetscCall(PetscSFDestroy(&overlapSF));
2392:   }

2394:   for (d = depth; d >= 0; d--) PetscCall(PlexLayerDestroy(&layers[d]));
2395:   PetscCall(PetscFree(layers));
2396:   PetscCall(ISDestroy(&strataPermutation));
2397:   PetscFunctionReturn(PETSC_SUCCESS);
2398: }

2400: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC)
2401: {
2402:   DMPlexStorageVersion version;
2403:   const char          *topologydm_name;
2404:   char                 group[PETSC_MAX_PATH_LEN];
2405:   PetscSF              sfwork = NULL;

2407:   PetscFunctionBegin;
2408:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2409:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2410:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
2411:     PetscCall(PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name));
2412:   } else {
2413:     PetscCall(PetscStrncpy(group, "/", sizeof(group)));
2414:   }
2415:   PetscCall(PetscViewerHDF5PushGroup(viewer, group));

2417:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topology"));
2418:   if (version->major < 3) {
2419:     PetscCall(DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork));
2420:   } else {
2421:     /* since DMPlexStorageVersion 3.0.0 */
2422:     PetscCall(DMPlexTopologyLoad_HDF5_Private(dm, viewer, &sfwork));
2423:   }
2424:   PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "topology" */

2426:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
2427:     DM          distdm;
2428:     PetscSF     distsf;
2429:     const char *distribution_name;
2430:     PetscBool   exists;

2432:     PetscCall(DMPlexDistributionGetName(dm, &distribution_name));
2433:     /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */
2434:     PetscCall(PetscViewerHDF5PushGroup(viewer, "distributions"));
2435:     PetscCall(PetscViewerHDF5HasGroup(viewer, NULL, &exists));
2436:     if (exists) {
2437:       PetscCall(PetscViewerHDF5PushGroup(viewer, distribution_name));
2438:       PetscCall(DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm));
2439:       if (distdm) {
2440:         PetscCall(DMPlexReplace_Internal(dm, &distdm));
2441:         PetscCall(PetscSFDestroy(&sfwork));
2442:         sfwork = distsf;
2443:       }
2444:       PetscCall(PetscViewerHDF5PopGroup(viewer)); /* distribution_name */
2445:     }
2446:     PetscCall(PetscViewerHDF5PopGroup(viewer)); /* "distributions" */
2447:   }
2448:   if (sfXC) {
2449:     *sfXC = sfwork;
2450:   } else {
2451:     PetscCall(PetscSFDestroy(&sfwork));
2452:   }

2454:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2455:   PetscFunctionReturn(PETSC_SUCCESS);
2456: }

2458: /* If the file is old, it not only has different path to the coordinates, but   */
2459: /* does not contain coordinateDMs, so must fall back to the old implementation. */
2460: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
2461: {
2462:   PetscSection coordSection;
2463:   Vec          coordinates;
2464:   PetscReal    lengthScale;
2465:   PetscInt     spatialDim, N, numVertices, vStart, vEnd, v;
2466:   PetscMPIInt  rank;

2468:   PetscFunctionBegin;
2469:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
2470:   /* Read geometry */
2471:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/geometry"));
2472:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), &coordinates));
2473:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "vertices"));
2474:   {
2475:     /* Force serial load */
2476:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N));
2477:     PetscCall(VecSetSizes(coordinates, rank == 0 ? N : 0, N));
2478:     PetscCall(VecSetBlockSize(coordinates, spatialDim));
2479:   }
2480:   PetscCall(VecLoad(coordinates, viewer));
2481:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2482:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
2483:   PetscCall(VecScale(coordinates, 1.0 / lengthScale));
2484:   PetscCall(VecGetLocalSize(coordinates, &numVertices));
2485:   PetscCall(VecGetBlockSize(coordinates, &spatialDim));
2486:   numVertices /= spatialDim;
2487:   /* Create coordinates */
2488:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
2489:   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);
2490:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2491:   PetscCall(PetscSectionSetNumFields(coordSection, 1));
2492:   PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spatialDim));
2493:   PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd));
2494:   for (v = vStart; v < vEnd; ++v) {
2495:     PetscCall(PetscSectionSetDof(coordSection, v, spatialDim));
2496:     PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spatialDim));
2497:   }
2498:   PetscCall(PetscSectionSetUp(coordSection));
2499:   PetscCall(DMSetCoordinates(dm, coordinates));
2500:   PetscCall(VecDestroy(&coordinates));
2501:   PetscFunctionReturn(PETSC_SUCCESS);
2502: }

2504: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
2505: {
2506:   DMPlexStorageVersion version;
2507:   DM                   cdm;
2508:   Vec                  coords;
2509:   PetscInt             blockSize;
2510:   PetscReal            lengthScale;
2511:   PetscSF              lsf;
2512:   const char          *topologydm_name;
2513:   char                *coordinatedm_name, *coordinates_name;

2515:   PetscFunctionBegin;
2516:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2517:   if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
2518:     PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
2519:     PetscFunctionReturn(PETSC_SUCCESS);
2520:   }
2521:   /* else: since DMPlexStorageVersion 2.0.0 */
2522:   PetscCheck(sfXC, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "PetscSF must be given for parallel load");
2523:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2524:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2525:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2526:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name));
2527:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name));
2528:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2529:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2530:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2531:   PetscCall(PetscObjectSetName((PetscObject)cdm, coordinatedm_name));
2532:   PetscCall(PetscFree(coordinatedm_name));
2533:   /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
2534:   PetscCall(DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf));
2535:   PetscCall(DMCreateLocalVector(cdm, &coords));
2536:   PetscCall(PetscObjectSetName((PetscObject)coords, coordinates_name));
2537:   PetscCall(PetscFree(coordinates_name));
2538:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE));
2539:   PetscCall(DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords));
2540:   PetscCall(PetscViewerPopFormat(viewer));
2541:   PetscCall(DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale));
2542:   PetscCall(VecScale(coords, 1.0 / lengthScale));
2543:   PetscCall(DMSetCoordinatesLocal(dm, coords));
2544:   PetscCall(VecGetBlockSize(coords, &blockSize));
2545:   PetscCall(DMSetCoordinateDim(dm, blockSize));
2546:   PetscCall(VecDestroy(&coords));
2547:   PetscCall(PetscSFDestroy(&lsf));
2548:   PetscFunctionReturn(PETSC_SUCCESS);
2549: }

2551: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
2552: {
2553:   DMPlexStorageVersion version;

2555:   PetscFunctionBegin;
2556:   PetscCall(PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version));
2557:   if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
2558:     PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL));
2559:     PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL));
2560:     PetscCall(DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer));
2561:   } else {
2562:     PetscSF sfXC;

2564:     /* since DMPlexStorageVersion 2.0.0 */
2565:     PetscCall(DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC));
2566:     PetscCall(DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC));
2567:     PetscCall(DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC));
2568:     PetscCall(PetscSFDestroy(&sfXC));
2569:   }
2570:   PetscFunctionReturn(PETSC_SUCCESS);
2571: }

2573: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
2574: {
2575:   MPI_Comm  comm;
2576:   PetscInt  pStart, pEnd, p, m;
2577:   PetscInt *goffs, *ilocal;
2578:   PetscBool rootIncludeConstraints, leafIncludeConstraints;

2580:   PetscFunctionBegin;
2581:   PetscCall(PetscObjectGetComm((PetscObject)leafSection, &comm));
2582:   PetscCall(PetscSectionGetChart(leafSection, &pStart, &pEnd));
2583:   PetscCall(PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints));
2584:   PetscCall(PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints));
2585:   if (rootIncludeConstraints && leafIncludeConstraints) PetscCall(PetscSectionGetStorageSize(leafSection, &m));
2586:   else PetscCall(PetscSectionGetConstrainedStorageSize(leafSection, &m));
2587:   PetscCall(PetscMalloc1(m, &ilocal));
2588:   PetscCall(PetscMalloc1(m, &goffs));
2589:   /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
2590:   /* for the top-level section (not for each field), so one must have   */
2591:   /* rootSection->pointMajor == PETSC_TRUE.                             */
2592:   PetscCheck(rootSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
2593:   /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
2594:   PetscCheck(leafSection->pointMajor, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for field major ordering");
2595:   for (p = pStart, m = 0; p < pEnd; ++p) {
2596:     PetscInt        dof, cdof, i, j, off, goff;
2597:     const PetscInt *cinds;

2599:     PetscCall(PetscSectionGetDof(leafSection, p, &dof));
2600:     if (dof < 0) continue;
2601:     goff = globalOffsets[p - pStart];
2602:     PetscCall(PetscSectionGetOffset(leafSection, p, &off));
2603:     PetscCall(PetscSectionGetConstraintDof(leafSection, p, &cdof));
2604:     PetscCall(PetscSectionGetConstraintIndices(leafSection, p, &cinds));
2605:     for (i = 0, j = 0; i < dof; ++i) {
2606:       PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]);

2608:       if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {
2609:         ilocal[m]  = off++;
2610:         goffs[m++] = goff++;
2611:       } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
2612:       else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
2613:       if (constrained) ++j;
2614:     }
2615:   }
2616:   PetscCall(PetscSFCreate(comm, sectionSF));
2617:   PetscCall(PetscSFSetFromOptions(*sectionSF));
2618:   PetscCall(PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs));
2619:   PetscCall(PetscFree(goffs));
2620:   PetscFunctionReturn(PETSC_SUCCESS);
2621: }

2623: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
2624: {
2625:   MPI_Comm     comm;
2626:   PetscMPIInt  size, rank;
2627:   const char  *topologydm_name;
2628:   const char  *sectiondm_name;
2629:   PetscSection sectionA, sectionB;
2630:   PetscInt     nX, n, i;
2631:   PetscSF      sfAB;

2633:   PetscFunctionBegin;
2634:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2635:   PetscCallMPI(MPI_Comm_size(comm, &size));
2636:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
2637:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2638:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
2639:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2640:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2641:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
2642:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
2643:   /* A: on-disk points                        */
2644:   /* X: list of global point numbers, [0, NX) */
2645:   /* B: plex points                           */
2646:   /* Load raw section (sectionA)              */
2647:   PetscCall(PetscSectionCreate(comm, &sectionA));
2648:   PetscCall(PetscSectionLoad(sectionA, viewer));
2649:   PetscCall(PetscSectionGetChart(sectionA, NULL, &n));
2650:   /* Create sfAB: A -> B */
2651:   #if defined(PETSC_USE_DEBUG)
2652:   {
2653:     PetscInt N, N1;

2655:     PetscCall(PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1));
2656:     PetscCall(MPIU_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm));
2657:     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);
2658:   }
2659:   #endif
2660:   {
2661:     IS              orderIS;
2662:     const PetscInt *gpoints;
2663:     PetscSF         sfXA, sfAX;
2664:     PetscLayout     layout;
2665:     PetscSFNode    *owners, *buffer;
2666:     PetscInt        nleaves;
2667:     PetscInt       *ilocal;
2668:     PetscSFNode    *iremote;

2670:     /* Create sfAX: A -> X */
2671:     PetscCall(ISCreate(comm, &orderIS));
2672:     PetscCall(PetscObjectSetName((PetscObject)orderIS, "order"));
2673:     PetscCall(PetscLayoutSetLocalSize(orderIS->map, n));
2674:     PetscCall(ISLoad(orderIS, viewer));
2675:     PetscCall(PetscLayoutCreate(comm, &layout));
2676:     PetscCall(PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL));
2677:     PetscCall(PetscLayoutSetLocalSize(layout, nX));
2678:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
2679:     PetscCall(PetscLayoutSetUp(layout));
2680:     PetscCall(PetscSFCreate(comm, &sfXA));
2681:     PetscCall(ISGetIndices(orderIS, &gpoints));
2682:     PetscCall(PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints));
2683:     PetscCall(ISRestoreIndices(orderIS, &gpoints));
2684:     PetscCall(ISDestroy(&orderIS));
2685:     PetscCall(PetscLayoutDestroy(&layout));
2686:     PetscCall(PetscMalloc1(n, &owners));
2687:     PetscCall(PetscMalloc1(nX, &buffer));
2688:     for (i = 0; i < n; ++i) {
2689:       owners[i].rank  = rank;
2690:       owners[i].index = i;
2691:     }
2692:     for (i = 0; i < nX; ++i) {
2693:       buffer[i].rank  = -1;
2694:       buffer[i].index = -1;
2695:     }
2696:     PetscCall(PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC));
2697:     PetscCall(PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC));
2698:     PetscCall(PetscSFDestroy(&sfXA));
2699:     PetscCall(PetscFree(owners));
2700:     for (i = 0, nleaves = 0; i < nX; ++i)
2701:       if (buffer[i].rank >= 0) nleaves++;
2702:     PetscCall(PetscMalloc1(nleaves, &ilocal));
2703:     PetscCall(PetscMalloc1(nleaves, &iremote));
2704:     for (i = 0, nleaves = 0; i < nX; ++i) {
2705:       if (buffer[i].rank >= 0) {
2706:         ilocal[nleaves]        = i;
2707:         iremote[nleaves].rank  = buffer[i].rank;
2708:         iremote[nleaves].index = buffer[i].index;
2709:         nleaves++;
2710:       }
2711:     }
2712:     PetscCall(PetscFree(buffer));
2713:     PetscCall(PetscSFCreate(comm, &sfAX));
2714:     PetscCall(PetscSFSetFromOptions(sfAX));
2715:     PetscCall(PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2716:     PetscCall(PetscSFCompose(sfAX, sfXB, &sfAB));
2717:     PetscCall(PetscSFDestroy(&sfAX));
2718:   }
2719:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2720:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2721:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2722:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2723:   /* Create plex section (sectionB) */
2724:   PetscCall(DMGetLocalSection(sectiondm, &sectionB));
2725:   if (lsf || gsf) {
2726:     PetscLayout layout;
2727:     PetscInt    M, m;
2728:     PetscInt   *offsetsA;
2729:     PetscBool   includesConstraintsA;

2731:     PetscCall(PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB));
2732:     PetscCall(PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA));
2733:     if (includesConstraintsA) PetscCall(PetscSectionGetStorageSize(sectionA, &m));
2734:     else PetscCall(PetscSectionGetConstrainedStorageSize(sectionA, &m));
2735:     PetscCall(MPIU_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm));
2736:     PetscCall(PetscLayoutCreate(comm, &layout));
2737:     PetscCall(PetscLayoutSetSize(layout, M));
2738:     PetscCall(PetscLayoutSetUp(layout));
2739:     if (lsf) {
2740:       PetscSF lsfABdata;

2742:       PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata));
2743:       *lsf = lsfABdata;
2744:     }
2745:     if (gsf) {
2746:       PetscSection gsectionB, gsectionB1;
2747:       PetscBool    includesConstraintsB;
2748:       PetscSF      gsfABdata, pointsf;

2750:       PetscCall(DMGetGlobalSection(sectiondm, &gsectionB1));
2751:       PetscCall(PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB));
2752:       PetscCall(DMGetPointSF(sectiondm, &pointsf));
2753:       PetscCall(PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB));
2754:       PetscCall(DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata));
2755:       PetscCall(PetscSectionDestroy(&gsectionB));
2756:       *gsf = gsfABdata;
2757:     }
2758:     PetscCall(PetscLayoutDestroy(&layout));
2759:     PetscCall(PetscFree(offsetsA));
2760:   } else {
2761:     PetscCall(PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB));
2762:   }
2763:   PetscCall(PetscSFDestroy(&sfAB));
2764:   PetscCall(PetscSectionDestroy(&sectionA));
2765:   PetscFunctionReturn(PETSC_SUCCESS);
2766: }

2768: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
2769: {
2770:   MPI_Comm           comm;
2771:   const char        *topologydm_name;
2772:   const char        *sectiondm_name;
2773:   const char        *vec_name;
2774:   Vec                vecA;
2775:   PetscInt           mA, m, bs;
2776:   const PetscInt    *ilocal;
2777:   const PetscScalar *src;
2778:   PetscScalar       *dest;

2780:   PetscFunctionBegin;
2781:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2782:   PetscCall(DMPlexGetHDF5Name_Private(dm, &topologydm_name));
2783:   PetscCall(PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name));
2784:   PetscCall(PetscObjectGetName((PetscObject)vec, &vec_name));
2785:   PetscCall(PetscViewerHDF5PushGroup(viewer, "topologies"));
2786:   PetscCall(PetscViewerHDF5PushGroup(viewer, topologydm_name));
2787:   PetscCall(PetscViewerHDF5PushGroup(viewer, "dms"));
2788:   PetscCall(PetscViewerHDF5PushGroup(viewer, sectiondm_name));
2789:   PetscCall(PetscViewerHDF5PushGroup(viewer, "vecs"));
2790:   PetscCall(PetscViewerHDF5PushGroup(viewer, vec_name));
2791:   PetscCall(VecCreate(comm, &vecA));
2792:   PetscCall(PetscObjectSetName((PetscObject)vecA, vec_name));
2793:   PetscCall(PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL));
2794:   /* Check consistency */
2795:   {
2796:     PetscSF  pointsf, pointsf1;
2797:     PetscInt m1, i, j;

2799:     PetscCall(DMGetPointSF(dm, &pointsf));
2800:     PetscCall(DMGetPointSF(sectiondm, &pointsf1));
2801:     PetscCheck(pointsf1 == pointsf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
2802:   #if defined(PETSC_USE_DEBUG)
2803:     {
2804:       PetscInt MA, MA1;

2806:       PetscCall(MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm));
2807:       PetscCall(PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1));
2808:       PetscCheck(MA1 == MA, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%" PetscInt_FMT ") != On-disk vector data size (%" PetscInt_FMT ")", MA, MA1);
2809:     }
2810:   #endif
2811:     PetscCall(VecGetLocalSize(vec, &m1));
2812:     PetscCheck(m1 >= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%" PetscInt_FMT ") < SF leaf size (%" PetscInt_FMT ")", m1, m);
2813:     for (i = 0; i < m; ++i) {
2814:       j = ilocal ? ilocal[i] : i;
2815:       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);
2816:     }
2817:   }
2818:   PetscCall(VecSetSizes(vecA, mA, PETSC_DECIDE));
2819:   PetscCall(VecLoad(vecA, viewer));
2820:   PetscCall(VecGetArrayRead(vecA, &src));
2821:   PetscCall(VecGetArray(vec, &dest));
2822:   PetscCall(PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
2823:   PetscCall(PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE));
2824:   PetscCall(VecRestoreArray(vec, &dest));
2825:   PetscCall(VecRestoreArrayRead(vecA, &src));
2826:   PetscCall(VecDestroy(&vecA));
2827:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs));
2828:   PetscCall(VecSetBlockSize(vec, bs));
2829:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2830:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2831:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2832:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2833:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2834:   PetscCall(PetscViewerHDF5PopGroup(viewer));
2835:   PetscFunctionReturn(PETSC_SUCCESS);
2836: }
2837: #endif