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

 13: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

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

 23:   PetscTokenCreate(str, '.', &t);
 24:   for (i = 0; i < 3; i++) {
 25:     PetscTokenFind(t, &ts);
 27:     PetscOptionsStringToInt(ts, &ti[i]);
 28:   }
 29:   PetscTokenFind(t, &ts);
 31:   PetscTokenDestroy(&t);
 32:   PetscNew(&v);
 33:   v->major    = ti[0];
 34:   v->minor    = ti[1];
 35:   v->subminor = ti[2];
 36:   PetscViewerCheckVersion_Private(viewer, v);
 37:   *version = v;
 38:   return 0;
 39: }

 41: static PetscErrorCode PetscViewerAttachVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion v)
 42: {
 43:   PetscContainer cont;

 45:   PetscContainerCreate(PetscObjectComm((PetscObject)viewer), &cont);
 46:   PetscContainerSetPointer(cont, v);
 47:   PetscContainerSetUserDestroy(cont, PetscContainerUserDestroyDefault);
 48:   PetscObjectCompose((PetscObject)viewer, key, (PetscObject)cont);
 49:   PetscContainerDestroy(&cont);
 50:   return 0;
 51: }

 53: static PetscErrorCode PetscViewerGetAttachedVersion_Private(PetscViewer viewer, const char key[], DMPlexStorageVersion *v)
 54: {
 55:   PetscContainer cont;

 57:   PetscObjectQuery((PetscObject)viewer, key, (PetscObject *)&cont);
 58:   *v = NULL;
 59:   if (cont) PetscContainerGetPointer(cont, (void **)v);
 60:   return 0;
 61: }

 63: /*
 64:   Version log:
 65:   1.0.0 legacy version (default if no "dmplex_storage_version" attribute found in file)
 66:   2.0.0 introduce versioning and multiple topologies
 67:   2.1.0 introduce distributions
 68: */
 69: static PetscErrorCode PetscViewerCheckVersion_Private(PetscViewer viewer, DMPlexStorageVersion version)
 70: {
 71:   PetscBool valid = PETSC_FALSE;

 73:   switch (version->major) {
 74:   case 1:
 75:     switch (version->minor) {
 76:     case 0:
 77:       switch (version->subminor) {
 78:       case 0:
 79:         valid = PETSC_TRUE;
 80:         break;
 81:       };
 82:       break;
 83:     };
 84:     break;
 85:   case 2:
 86:     switch (version->minor) {
 87:     case 0:
 88:       switch (version->subminor) {
 89:       case 0:
 90:         valid = PETSC_TRUE;
 91:         break;
 92:       };
 93:       break;
 94:     case 1:
 95:       switch (version->subminor) {
 96:       case 0:
 97:         valid = PETSC_TRUE;
 98:         break;
 99:       };
100:       break;
101:     };
102:     break;
103:   }
105:   return 0;
106: }

108: static inline PetscBool DMPlexStorageVersionGE(DMPlexStorageVersion version, int major, int minor, int subminor)
109: {
110:   return (PetscBool)((version->major == major && version->minor == minor && version->subminor >= subminor) || (version->major == major && version->minor > minor) || (version->major > major));
111: }

113: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer viewer, DMPlexStorageVersion *version)
114: {
115:   const char ATTR_NAME[] = "dmplex_storage_version";
116:   PetscBool  fileHasVersion;
117:   char       optVersion[16], fileVersion[16];

121:   PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, version);
122:   if (*version) return 0;

124:   PetscStrcpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE);
125:   PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion);
126:   if (fileHasVersion) {
127:     char *tmp;

129:     PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, NULL, &tmp);
130:     PetscStrcpy(fileVersion, tmp);
131:     PetscFree(tmp);
132:   }
133:   PetscStrcpy(optVersion, fileVersion);
134:   PetscObjectOptionsBegin((PetscObject)viewer);
135:   PetscOptionsString("-dm_plex_view_hdf5_storage_version", "DMPlex HDF5 viewer storage version", NULL, optVersion, optVersion, sizeof(optVersion), NULL);
136:   PetscOptionsEnd();
137:   if (!fileHasVersion) {
138:     PetscViewerHDF5WriteAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, optVersion);
139:   } else {
140:     PetscBool flg;

142:     PetscStrcmp(fileVersion, optVersion, &flg);
144:   }
145:   PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);
146:   PetscViewerParseVersion_Private(viewer, optVersion, version);
147:   PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_WRITING_KEY, *version);
148:   return 0;
149: }

151: PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer viewer, DMPlexStorageVersion *version)
152: {
153:   const char ATTR_NAME[] = "dmplex_storage_version";
154:   char      *defaultVersion;
155:   char      *versionString;

159:   PetscViewerGetAttachedVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, version);
160:   if (*version) return 0;

162:   //TODO string HDF5 attribute handling is terrible and should be redesigned
163:   PetscStrallocpy("1.0.0", &defaultVersion);
164:   PetscViewerHDF5ReadAttribute(viewer, "/", ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString);
165:   PetscViewerParseVersion_Private(viewer, versionString, version);
166:   PetscViewerAttachVersion_Private(viewer, DMPLEX_STORAGE_VERSION_READING_KEY, *version);
167:   PetscFree(versionString);
168:   PetscFree(defaultVersion);
169:   return 0;
170: }

172: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
173: {
174:   if (((PetscObject)dm)->name) {
175:     PetscObjectGetName((PetscObject)dm, name);
176:   } else {
177:     *name = "plex";
178:   }
179:   return 0;
180: }

182: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
183: {
184:   Vec         stamp;
185:   PetscMPIInt rank;

187:   if (seqnum < 0) return 0;
188:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
189:   VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp);
190:   VecSetBlockSize(stamp, 1);
191:   PetscObjectSetName((PetscObject)stamp, seqname);
192:   if (rank == 0) {
193:     PetscReal timeScale;
194:     PetscBool istime;

196:     PetscStrncmp(seqname, "time", 5, &istime);
197:     if (istime) {
198:       DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale);
199:       value *= timeScale;
200:     }
201:     VecSetValue(stamp, 0, value, INSERT_VALUES);
202:   }
203:   VecAssemblyBegin(stamp);
204:   VecAssemblyEnd(stamp);
205:   PetscViewerHDF5PushGroup(viewer, "/");
206:   PetscViewerHDF5PushTimestepping(viewer);
207:   PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
208:   VecView(stamp, viewer);
209:   PetscViewerHDF5PopTimestepping(viewer);
210:   PetscViewerHDF5PopGroup(viewer);
211:   VecDestroy(&stamp);
212:   return 0;
213: }

215: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
216: {
217:   Vec         stamp;
218:   PetscMPIInt rank;

220:   if (seqnum < 0) return 0;
221:   MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank);
222:   VecCreateMPI(PetscObjectComm((PetscObject)viewer), rank ? 0 : 1, 1, &stamp);
223:   VecSetBlockSize(stamp, 1);
224:   PetscObjectSetName((PetscObject)stamp, seqname);
225:   PetscViewerHDF5PushGroup(viewer, "/");
226:   PetscViewerHDF5PushTimestepping(viewer);
227:   PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
228:   VecLoad(stamp, viewer);
229:   PetscViewerHDF5PopTimestepping(viewer);
230:   PetscViewerHDF5PopGroup(viewer);
231:   if (rank == 0) {
232:     const PetscScalar *a;
233:     PetscReal          timeScale;
234:     PetscBool          istime;

236:     VecGetArrayRead(stamp, &a);
237:     *value = a[0];
238:     VecRestoreArrayRead(stamp, &a);
239:     PetscStrncmp(seqname, "time", 5, &istime);
240:     if (istime) {
241:       DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale);
242:       *value /= timeScale;
243:     }
244:   }
245:   VecDestroy(&stamp);
246:   return 0;
247: }

249: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
250: {
251:   IS              cutcells = NULL;
252:   const PetscInt *cutc;
253:   PetscInt        cellHeight, vStart, vEnd, cStart, cEnd, c;

255:   if (!cutLabel) return 0;
256:   DMPlexGetVTKCellHeight(dm, &cellHeight);
257:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
258:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
259:   /* Label vertices that should be duplicated */
260:   DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
261:   DMLabelGetStratumIS(cutLabel, 2, &cutcells);
262:   if (cutcells) {
263:     PetscInt n;

265:     ISGetIndices(cutcells, &cutc);
266:     ISGetLocalSize(cutcells, &n);
267:     for (c = 0; c < n; ++c) {
268:       if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
269:         PetscInt *closure = NULL;
270:         PetscInt  closureSize, cl, value;

272:         DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
273:         for (cl = 0; cl < closureSize * 2; cl += 2) {
274:           if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
275:             DMLabelGetValue(cutLabel, closure[cl], &value);
276:             if (value == 1) DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
277:           }
278:         }
279:         DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
280:       }
281:     }
282:     ISRestoreIndices(cutcells, &cutc);
283:   }
284:   ISDestroy(&cutcells);
285:   return 0;
286: }

288: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
289: {
290:   DM                      dm;
291:   DM                      dmBC;
292:   PetscSection            section, sectionGlobal;
293:   Vec                     gv;
294:   const char             *name;
295:   PetscViewerVTKFieldType ft;
296:   PetscViewerFormat       format;
297:   PetscInt                seqnum;
298:   PetscReal               seqval;
299:   PetscBool               isseq;

301:   PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
302:   VecGetDM(v, &dm);
303:   DMGetLocalSection(dm, &section);
304:   DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
305:   DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar)seqval, viewer);
306:   if (seqnum >= 0) {
307:     PetscViewerHDF5PushTimestepping(viewer);
308:     PetscViewerHDF5SetTimestep(viewer, seqnum);
309:   }
310:   PetscViewerGetFormat(viewer, &format);
311:   DMGetOutputDM(dm, &dmBC);
312:   DMGetGlobalSection(dmBC, &sectionGlobal);
313:   DMGetGlobalVector(dmBC, &gv);
314:   PetscObjectGetName((PetscObject)v, &name);
315:   PetscObjectSetName((PetscObject)gv, name);
316:   DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
317:   DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
318:   PetscObjectTypeCompare((PetscObject)gv, VECSEQ, &isseq);
319:   if (isseq) VecView_Seq(gv, viewer);
320:   else VecView_MPI(gv, viewer);
321:   if (format == PETSC_VIEWER_HDF5_VIZ) {
322:     /* Output visualization representation */
323:     PetscInt numFields, f;
324:     DMLabel  cutLabel, cutVertexLabel = NULL;

326:     PetscSectionGetNumFields(section, &numFields);
327:     DMGetLabel(dm, "periodic_cut", &cutLabel);
328:     for (f = 0; f < numFields; ++f) {
329:       Vec         subv;
330:       IS          is;
331:       const char *fname, *fgroup, *componentName;
332:       char        subname[PETSC_MAX_PATH_LEN];
333:       PetscInt    pStart, pEnd, Nc, c;

335:       DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
336:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
337:       PetscSectionGetFieldName(section, f, &fname);
338:       if (!fname || ft == PETSC_VTK_INVALID) continue;
339:       PetscViewerHDF5PushGroup(viewer, fgroup);
340:       if (cutLabel) {
341:         const PetscScalar *ga;
342:         PetscScalar       *suba;
343:         PetscInt           gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;

345:         DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
346:         PetscSectionGetFieldComponents(section, f, &Nc);
347:         for (p = pStart; p < pEnd; ++p) {
348:           PetscInt gdof, fdof = 0, val;

350:           PetscSectionGetDof(sectionGlobal, p, &gdof);
351:           if (gdof > 0) PetscSectionGetFieldDof(section, p, f, &fdof);
352:           subSize += fdof;
353:           DMLabelGetValue(cutVertexLabel, p, &val);
354:           if (val == 1) extSize += fdof;
355:         }
356:         VecCreate(PetscObjectComm((PetscObject)gv), &subv);
357:         VecSetSizes(subv, subSize + extSize, PETSC_DETERMINE);
358:         VecSetBlockSize(subv, Nc);
359:         VecSetType(subv, VECSTANDARD);
360:         VecGetOwnershipRange(gv, &gstart, NULL);
361:         VecGetArrayRead(gv, &ga);
362:         VecGetArray(subv, &suba);
363:         for (p = pStart; p < pEnd; ++p) {
364:           PetscInt gdof, goff, val;

366:           PetscSectionGetDof(sectionGlobal, p, &gdof);
367:           if (gdof > 0) {
368:             PetscInt fdof, fc, f2, poff = 0;

370:             PetscSectionGetOffset(sectionGlobal, p, &goff);
371:             /* Can get rid of this loop by storing field information in the global section */
372:             for (f2 = 0; f2 < f; ++f2) {
373:               PetscSectionGetFieldDof(section, p, f2, &fdof);
374:               poff += fdof;
375:             }
376:             PetscSectionGetFieldDof(section, p, f, &fdof);
377:             for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff + poff + fc - gstart];
378:             DMLabelGetValue(cutVertexLabel, p, &val);
379:             if (val == 1) {
380:               for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize + newOff] = ga[goff + poff + fc - gstart];
381:             }
382:           }
383:         }
384:         VecRestoreArrayRead(gv, &ga);
385:         VecRestoreArray(subv, &suba);
386:         DMLabelDestroy(&cutVertexLabel);
387:       } else {
388:         PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
389:       }
390:       PetscStrncpy(subname, name, sizeof(subname));
391:       PetscStrlcat(subname, "_", sizeof(subname));
392:       PetscStrlcat(subname, fname, sizeof(subname));
393:       PetscObjectSetName((PetscObject)subv, subname);
394:       if (isseq) VecView_Seq(subv, viewer);
395:       else VecView_MPI(subv, viewer);
396:       if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
397:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "vector");
398:       } else {
399:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, "vector_field_type", PETSC_STRING, "scalar");
400:       }

402:       /* Output the component names in the field if available */
403:       PetscSectionGetFieldComponents(section, f, &Nc);
404:       for (c = 0; c < Nc; ++c) {
405:         char componentNameLabel[PETSC_MAX_PATH_LEN];
406:         PetscSectionGetComponentName(section, f, c, &componentName);
407:         PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%" PetscInt_FMT, c);
408:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)subv, componentNameLabel, PETSC_STRING, componentName);
409:       }

411:       if (cutLabel) VecDestroy(&subv);
412:       else PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
413:       PetscViewerHDF5PopGroup(viewer);
414:     }
415:   }
416:   if (seqnum >= 0) PetscViewerHDF5PopTimestepping(viewer);
417:   DMRestoreGlobalVector(dmBC, &gv);
418:   return 0;
419: }

421: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
422: {
423:   DM          dm;
424:   Vec         locv;
425:   PetscObject isZero;
426:   const char *name;
427:   PetscReal   time;

429:   VecGetDM(v, &dm);
430:   DMGetLocalVector(dm, &locv);
431:   PetscObjectGetName((PetscObject)v, &name);
432:   PetscObjectSetName((PetscObject)locv, name);
433:   PetscObjectQuery((PetscObject)v, "__Vec_bc_zero__", &isZero);
434:   PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", isZero);
435:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
436:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
437:   DMGetOutputSequenceNumber(dm, NULL, &time);
438:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
439:   PetscViewerHDF5PushGroup(viewer, "/fields");
440:   PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
441:   VecView_Plex_Local_HDF5_Internal(locv, viewer);
442:   PetscViewerPopFormat(viewer);
443:   PetscViewerHDF5PopGroup(viewer);
444:   PetscObjectCompose((PetscObject)locv, "__Vec_bc_zero__", NULL);
445:   DMRestoreLocalVector(dm, &locv);
446:   return 0;
447: }

449: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
450: {
451:   PetscBool isseq;

453:   PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
454:   PetscViewerHDF5PushGroup(viewer, "/fields");
455:   if (isseq) VecView_Seq(v, viewer);
456:   else VecView_MPI(v, viewer);
457:   PetscViewerHDF5PopGroup(viewer);
458:   return 0;
459: }

461: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
462: {
463:   DM          dm;
464:   Vec         locv;
465:   const char *name;
466:   PetscInt    seqnum;

468:   VecGetDM(v, &dm);
469:   DMGetLocalVector(dm, &locv);
470:   PetscObjectGetName((PetscObject)v, &name);
471:   PetscObjectSetName((PetscObject)locv, name);
472:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
473:   PetscViewerHDF5PushGroup(viewer, "/fields");
474:   if (seqnum >= 0) {
475:     PetscViewerHDF5PushTimestepping(viewer);
476:     PetscViewerHDF5SetTimestep(viewer, seqnum);
477:   }
478:   VecLoad_Plex_Local(locv, viewer);
479:   if (seqnum >= 0) PetscViewerHDF5PopTimestepping(viewer);
480:   PetscViewerHDF5PopGroup(viewer);
481:   DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
482:   DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
483:   DMRestoreLocalVector(dm, &locv);
484:   return 0;
485: }

487: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
488: {
489:   DM       dm;
490:   PetscInt seqnum;

492:   VecGetDM(v, &dm);
493:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
494:   PetscViewerHDF5PushGroup(viewer, "/fields");
495:   if (seqnum >= 0) {
496:     PetscViewerHDF5PushTimestepping(viewer);
497:     PetscViewerHDF5SetTimestep(viewer, seqnum);
498:   }
499:   VecLoad_Default(v, viewer);
500:   if (seqnum >= 0) PetscViewerHDF5PopTimestepping(viewer);
501:   PetscViewerHDF5PopGroup(viewer);
502:   return 0;
503: }

505: static PetscErrorCode DMPlexDistributionView_HDF5_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
506: {
507:   MPI_Comm           comm;
508:   PetscMPIInt        size, rank;
509:   const char        *topologydm_name, *distribution_name;
510:   const PetscInt    *gpoint;
511:   PetscInt           pStart, pEnd, p;
512:   PetscSF            pointSF;
513:   PetscInt           nroots, nleaves;
514:   const PetscInt    *ilocal;
515:   const PetscSFNode *iremote;
516:   IS                 chartSizesIS, ownersIS, gpointsIS;
517:   PetscInt          *chartSize, *owners, *gpoints;

519:   PetscObjectGetComm((PetscObject)dm, &comm);
520:   MPI_Comm_size(comm, &size);
521:   MPI_Comm_rank(comm, &rank);
522:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
523:   DMPlexDistributionGetName(dm, &distribution_name);
524:   if (!distribution_name) return 0;
525:   PetscViewerHDF5WriteAttribute(viewer, NULL, "comm_size", PETSC_INT, (void *)&size);
526:   ISGetIndices(globalPointNumbers, &gpoint);
527:   DMPlexGetChart(dm, &pStart, &pEnd);
528:   PetscMalloc1(1, &chartSize);
529:   *chartSize = pEnd - pStart;
530:   PetscMalloc1(*chartSize, &owners);
531:   PetscMalloc1(*chartSize, &gpoints);
532:   DMGetPointSF(dm, &pointSF);
533:   PetscSFGetGraph(pointSF, &nroots, &nleaves, &ilocal, &iremote);
534:   for (p = pStart; p < pEnd; ++p) {
535:     PetscInt gp = gpoint[p - pStart];

537:     owners[p - pStart]  = rank;
538:     gpoints[p - pStart] = (gp < 0 ? -(gp + 1) : gp);
539:   }
540:   for (p = 0; p < nleaves; ++p) {
541:     PetscInt ilocalp = (ilocal ? ilocal[p] : p);

543:     owners[ilocalp] = iremote[p].rank;
544:   }
545:   ISCreateGeneral(comm, 1, chartSize, PETSC_OWN_POINTER, &chartSizesIS);
546:   ISCreateGeneral(comm, *chartSize, owners, PETSC_OWN_POINTER, &ownersIS);
547:   ISCreateGeneral(comm, *chartSize, gpoints, PETSC_OWN_POINTER, &gpointsIS);
548:   PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes");
549:   PetscObjectSetName((PetscObject)ownersIS, "owners");
550:   PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers");
551:   ISView(chartSizesIS, viewer);
552:   ISView(ownersIS, viewer);
553:   ISView(gpointsIS, viewer);
554:   ISDestroy(&chartSizesIS);
555:   ISDestroy(&ownersIS);
556:   ISDestroy(&gpointsIS);
557:   ISRestoreIndices(globalPointNumbers, &gpoint);
558:   return 0;
559: }

561: 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[])
562: {
563:   IS              coneSizesIS, conesIS, orientationsIS;
564:   PetscInt       *coneSizes, *cones, *orientations;
565:   const PetscInt *gpoint;
566:   PetscInt        nPoints = 0, conesSize = 0;
567:   PetscInt        p, c, s;
568:   MPI_Comm        comm;

570:   PetscObjectGetComm((PetscObject)dm, &comm);
571:   ISGetIndices(globalPointNumbers, &gpoint);
572:   for (p = pStart; p < pEnd; ++p) {
573:     if (gpoint[p] >= 0) {
574:       PetscInt coneSize;

576:       DMPlexGetConeSize(dm, p, &coneSize);
577:       nPoints += 1;
578:       conesSize += coneSize;
579:     }
580:   }
581:   PetscMalloc1(nPoints, &coneSizes);
582:   PetscMalloc1(conesSize, &cones);
583:   PetscMalloc1(conesSize, &orientations);
584:   for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
585:     if (gpoint[p] >= 0) {
586:       const PetscInt *cone, *ornt;
587:       PetscInt        coneSize, cp;

589:       DMPlexGetConeSize(dm, p, &coneSize);
590:       DMPlexGetCone(dm, p, &cone);
591:       DMPlexGetConeOrientation(dm, p, &ornt);
592:       coneSizes[s] = coneSize;
593:       for (cp = 0; cp < coneSize; ++cp, ++c) {
594:         cones[c]        = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]] + 1) : gpoint[cone[cp]];
595:         orientations[c] = ornt[cp];
596:       }
597:       ++s;
598:     }
599:   }
602:   ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS);
603:   ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS);
604:   ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS);
605:   PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName);
606:   PetscObjectSetName((PetscObject)conesIS, conesName);
607:   PetscObjectSetName((PetscObject)orientationsIS, orientationsName);
608:   ISView(coneSizesIS, viewer);
609:   ISView(conesIS, viewer);
610:   ISView(orientationsIS, viewer);
611:   ISDestroy(&coneSizesIS);
612:   ISDestroy(&conesIS);
613:   ISDestroy(&orientationsIS);
614:   if (pointsName) {
615:     IS        pointsIS;
616:     PetscInt *points;

618:     PetscMalloc1(nPoints, &points);
619:     for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
620:       if (gpoint[p] >= 0) {
621:         points[s] = gpoint[p];
622:         ++s;
623:       }
624:     }
625:     ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS);
626:     PetscObjectSetName((PetscObject)pointsIS, pointsName);
627:     ISView(pointsIS, viewer);
628:     ISDestroy(&pointsIS);
629:   }
630:   ISRestoreIndices(globalPointNumbers, &gpoint);
631:   return 0;
632: }

634: static PetscErrorCode DMPlexTopologyView_HDF5_Legacy_Private(DM dm, IS globalPointNumbers, PetscViewer viewer)
635: {
636:   const char *pointsName, *coneSizesName, *conesName, *orientationsName;
637:   PetscInt    pStart, pEnd, dim;

639:   pointsName       = "order";
640:   coneSizesName    = "cones";
641:   conesName        = "cells";
642:   orientationsName = "orientation";
643:   DMPlexGetChart(dm, &pStart, &pEnd);
644:   DMPlexTopologyView_HDF5_Inner_Private(dm, globalPointNumbers, viewer, pStart, pEnd, pointsName, coneSizesName, conesName, orientationsName);
645:   DMGetDimension(dm, &dim);
646:   PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *)&dim);
647:   return 0;
648: }

650: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
651: {
652:   DMPlexStorageVersion version;
653:   const char          *topologydm_name;
654:   char                 group[PETSC_MAX_PATH_LEN];

656:   PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version);
657:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
658:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
659:     PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name);
660:   } else {
661:     PetscStrcpy(group, "/");
662:   }
663:   PetscViewerHDF5PushGroup(viewer, group);

665:   PetscViewerHDF5PushGroup(viewer, "topology");
666:   DMPlexTopologyView_HDF5_Legacy_Private(dm, globalPointNumbers, viewer);
667:   PetscViewerHDF5PopGroup(viewer); /* "topology" */

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

672:     DMPlexDistributionGetName(dm, &distribution_name);
673:     PetscViewerHDF5PushGroup(viewer, "distributions");
674:     PetscViewerHDF5WriteGroup(viewer, NULL);
675:     PetscViewerHDF5PushGroup(viewer, distribution_name);
676:     DMPlexDistributionView_HDF5_Private(dm, globalPointNumbers, viewer);
677:     PetscViewerHDF5PopGroup(viewer); /* distribution_name */
678:     PetscViewerHDF5PopGroup(viewer); /* "distributions" */
679:   }

681:   PetscViewerHDF5PopGroup(viewer);
682:   return 0;
683: }

685: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
686: {
687:   PetscSF         sfPoint;
688:   DMLabel         cutLabel, cutVertexLabel         = NULL;
689:   IS              globalVertexNumbers, cutvertices = NULL;
690:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
691:   PetscInt       *vertices;
692:   PetscInt        conesSize = 0;
693:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;

695:   *numCorners = 0;
696:   DMGetDimension(dm, &dim);
697:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
698:   ISGetIndices(globalCellNumbers, &gcell);

700:   for (cell = cStart; cell < cEnd; ++cell) {
701:     PetscInt *closure = NULL;
702:     PetscInt  closureSize, v, Nc = 0;

704:     if (gcell[cell] < 0) continue;
705:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
706:     for (v = 0; v < closureSize * 2; v += 2) {
707:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
708:     }
709:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
710:     conesSize += Nc;
711:     if (!numCornersLocal) numCornersLocal = Nc;
712:     else if (numCornersLocal != Nc) numCornersLocal = 1;
713:   }
714:   MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
716:   /* Handle periodic cuts by identifying vertices which should be duplicated */
717:   DMGetLabel(dm, "periodic_cut", &cutLabel);
718:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
719:   if (cutVertexLabel) DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);
720:   if (cutvertices) {
721:     ISGetIndices(cutvertices, &cutverts);
722:     ISGetLocalSize(cutvertices, &vExtra);
723:   }
724:   DMGetPointSF(dm, &sfPoint);
725:   if (cutLabel) {
726:     const PetscInt    *ilocal;
727:     const PetscSFNode *iremote;
728:     PetscInt           nroots, nleaves;

730:     PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
731:     if (nleaves < 0) {
732:       PetscObjectReference((PetscObject)sfPoint);
733:     } else {
734:       PetscSFCreate(PetscObjectComm((PetscObject)sfPoint), &sfPoint);
735:       PetscSFSetGraph(sfPoint, nroots + vExtra, nleaves, (PetscInt *)ilocal, PETSC_COPY_VALUES, (PetscSFNode *)iremote, PETSC_COPY_VALUES);
736:     }
737:   } else {
738:     PetscObjectReference((PetscObject)sfPoint);
739:   }
740:   /* Number all vertices */
741:   DMPlexCreateNumbering_Plex(dm, vStart, vEnd + vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
742:   PetscSFDestroy(&sfPoint);
743:   /* Create cones */
744:   ISGetIndices(globalVertexNumbers, &gvertex);
745:   PetscMalloc1(conesSize, &vertices);
746:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
747:     PetscInt *closure = NULL;
748:     PetscInt  closureSize, Nc = 0, p, value = -1;
749:     PetscBool replace;

751:     if (gcell[cell] < 0) continue;
752:     if (cutLabel) DMLabelGetValue(cutLabel, cell, &value);
753:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
754:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
755:     for (p = 0; p < closureSize * 2; p += 2) {
756:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) closure[Nc++] = closure[p];
757:     }
758:     DMPlexReorderCell(dm, cell, closure);
759:     for (p = 0; p < Nc; ++p) {
760:       PetscInt nv, gv = gvertex[closure[p] - vStart];

762:       if (replace) {
763:         PetscFindInt(closure[p], vExtra, cutverts, &nv);
764:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
765:       }
766:       vertices[v++] = gv < 0 ? -(gv + 1) : gv;
767:     }
768:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
769:   }
770:   ISRestoreIndices(globalVertexNumbers, &gvertex);
771:   ISDestroy(&globalVertexNumbers);
772:   ISRestoreIndices(globalCellNumbers, &gcell);
773:   if (cutvertices) ISRestoreIndices(cutvertices, &cutverts);
774:   ISDestroy(&cutvertices);
775:   DMLabelDestroy(&cutVertexLabel);
777:   ISCreateGeneral(PetscObjectComm((PetscObject)dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
778:   PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
779:   PetscObjectSetName((PetscObject)*cellIS, "cells");
780:   return 0;
781: }

783: static PetscErrorCode DMPlexTopologyView_HDF5_XDMF_Private(DM dm, IS globalCellNumbers, PetscViewer viewer)
784: {
785:   DM       cdm;
786:   DMLabel  depthLabel, ctLabel;
787:   IS       cellIS;
788:   PetscInt dim, depth, cellHeight, c, n = 0;

790:   PetscViewerHDF5PushGroup(viewer, "/viz");
791:   PetscViewerHDF5WriteGroup(viewer, NULL);

793:   PetscViewerHDF5PopGroup(viewer);
794:   DMGetDimension(dm, &dim);
795:   DMPlexGetDepth(dm, &depth);
796:   DMGetCoordinateDM(dm, &cdm);
797:   DMPlexGetVTKCellHeight(dm, &cellHeight);
798:   DMPlexGetDepthLabel(dm, &depthLabel);
799:   DMPlexGetCellTypeLabel(dm, &ctLabel);
800:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
801:     const DMPolytopeType ict = (DMPolytopeType)c;
802:     PetscInt             pStart, pEnd, dep, numCorners;
803:     PetscBool            output = PETSC_FALSE, doOutput;

805:     if (ict == DM_POLYTOPE_FV_GHOST) continue;
806:     DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
807:     if (pStart >= 0) {
808:       DMLabelGetValue(depthLabel, pStart, &dep);
809:       if (dep == depth - cellHeight) output = PETSC_TRUE;
810:     }
811:     MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)dm));
812:     if (!doOutput) continue;
813:     CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners, &cellIS);
814:     if (!n) {
815:       PetscViewerHDF5PushGroup(viewer, "/viz/topology");
816:     } else {
817:       char group[PETSC_MAX_PATH_LEN];

819:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%" PetscInt_FMT, n);
820:       PetscViewerHDF5PushGroup(viewer, group);
821:     }
822:     ISView(cellIS, viewer);
823:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_corners", PETSC_INT, (void *)&numCorners);
824:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)cellIS, "cell_dim", PETSC_INT, (void *)&dim);
825:     ISDestroy(&cellIS);
826:     PetscViewerHDF5PopGroup(viewer);
827:     ++n;
828:   }
829:   return 0;
830: }

832: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
833: {
834:   DM        cdm;
835:   Vec       coordinates, newcoords;
836:   PetscReal lengthScale;
837:   PetscInt  m, M, bs;

839:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
840:   DMGetCoordinateDM(dm, &cdm);
841:   DMGetCoordinates(dm, &coordinates);
842:   VecCreate(PetscObjectComm((PetscObject)coordinates), &newcoords);
843:   PetscObjectSetName((PetscObject)newcoords, "vertices");
844:   VecGetSize(coordinates, &M);
845:   VecGetLocalSize(coordinates, &m);
846:   VecSetSizes(newcoords, m, M);
847:   VecGetBlockSize(coordinates, &bs);
848:   VecSetBlockSize(newcoords, bs);
849:   VecSetType(newcoords, VECSTANDARD);
850:   VecCopy(coordinates, newcoords);
851:   VecScale(newcoords, lengthScale);
852:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
853:   PetscViewerHDF5PushGroup(viewer, "/geometry");
854:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
855:   VecView(newcoords, viewer);
856:   PetscViewerPopFormat(viewer);
857:   PetscViewerHDF5PopGroup(viewer);
858:   VecDestroy(&newcoords);
859:   return 0;
860: }

862: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
863: {
864:   DM          cdm;
865:   Vec         coords, newcoords;
866:   PetscInt    m, M, bs;
867:   PetscReal   lengthScale;
868:   const char *topologydm_name, *coordinatedm_name, *coordinates_name;

870:   {
871:     PetscViewerFormat    format;
872:     DMPlexStorageVersion version;

874:     PetscViewerGetFormat(viewer, &format);
875:     PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version);
876:     if (!DMPlexStorageVersionGE(version, 2, 0, 0) || format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ) {
877:       DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer);
878:       return 0;
879:     }
880:   }
881:   /* since 2.0.0 */
882:   DMGetCoordinateDM(dm, &cdm);
883:   DMGetCoordinates(dm, &coords);
884:   PetscObjectGetName((PetscObject)cdm, &coordinatedm_name);
885:   PetscObjectGetName((PetscObject)coords, &coordinates_name);
886:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
887:   PetscViewerHDF5PushGroup(viewer, "topologies");
888:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
889:   PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name);
890:   PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name);
891:   PetscViewerHDF5PopGroup(viewer);
892:   PetscViewerHDF5PopGroup(viewer);
893:   DMPlexSectionView(dm, viewer, cdm);
894:   VecCreate(PetscObjectComm((PetscObject)coords), &newcoords);
895:   PetscObjectSetName((PetscObject)newcoords, coordinates_name);
896:   VecGetSize(coords, &M);
897:   VecGetLocalSize(coords, &m);
898:   VecSetSizes(newcoords, m, M);
899:   VecGetBlockSize(coords, &bs);
900:   VecSetBlockSize(newcoords, bs);
901:   VecSetType(newcoords, VECSTANDARD);
902:   VecCopy(coords, newcoords);
903:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
904:   VecScale(newcoords, lengthScale);
905:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
906:   DMPlexGlobalVectorView(dm, viewer, cdm, newcoords);
907:   PetscViewerPopFormat(viewer);
908:   VecDestroy(&newcoords);
909:   return 0;
910: }

912: static PetscErrorCode DMPlexCoordinatesView_HDF5_XDMF_Private(DM dm, PetscViewer viewer)
913: {
914:   DM               cdm;
915:   Vec              coordinatesLocal, newcoords;
916:   PetscSection     cSection, cGlobalSection;
917:   PetscScalar     *coords, *ncoords;
918:   DMLabel          cutLabel, cutVertexLabel = NULL;
919:   const PetscReal *L;
920:   PetscReal        lengthScale;
921:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
922:   PetscBool        localized, embedded;

924:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
925:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
926:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
927:   VecGetBlockSize(coordinatesLocal, &bs);
928:   DMGetCoordinatesLocalized(dm, &localized);
929:   if (localized == PETSC_FALSE) return 0;
930:   DMGetPeriodicity(dm, NULL, NULL, &L);
931:   DMGetCoordinateDM(dm, &cdm);
932:   DMGetLocalSection(cdm, &cSection);
933:   DMGetGlobalSection(cdm, &cGlobalSection);
934:   DMGetLabel(dm, "periodic_cut", &cutLabel);
935:   N = 0;

937:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
938:   VecCreate(PetscObjectComm((PetscObject)dm), &newcoords);
939:   PetscSectionGetDof(cSection, vStart, &dof);
940:   PetscPrintf(PETSC_COMM_SELF, "DOF: %" PetscInt_FMT "\n", dof);
941:   embedded = (PetscBool)(L && dof == 2 && !cutLabel);
942:   if (cutVertexLabel) {
943:     DMLabelGetStratumSize(cutVertexLabel, 1, &v);
944:     N += dof * v;
945:   }
946:   for (v = vStart; v < vEnd; ++v) {
947:     PetscSectionGetDof(cGlobalSection, v, &dof);
948:     if (dof < 0) continue;
949:     if (embedded) N += dof + 1;
950:     else N += dof;
951:   }
952:   if (embedded) VecSetBlockSize(newcoords, bs + 1);
953:   else VecSetBlockSize(newcoords, bs);
954:   VecSetSizes(newcoords, N, PETSC_DETERMINE);
955:   VecSetType(newcoords, VECSTANDARD);
956:   VecGetArray(coordinatesLocal, &coords);
957:   VecGetArray(newcoords, &ncoords);
958:   coordSize = 0;
959:   for (v = vStart; v < vEnd; ++v) {
960:     PetscSectionGetDof(cGlobalSection, v, &dof);
961:     PetscSectionGetOffset(cSection, v, &off);
962:     if (dof < 0) continue;
963:     if (embedded) {
964:       if (L && (L[0] > 0.0) && (L[1] > 0.0)) {
965:         PetscReal theta, phi, r, R;
966:         /* XY-periodic */
967:         /* Suppose its an y-z circle, then
968:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
969:            and the circle in that plane is
970:              \hat r cos(phi) + \hat x sin(phi) */
971:         theta                = 2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1];
972:         phi                  = 2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0];
973:         r                    = L[0] / (2.0 * PETSC_PI * 2.0 * L[1]);
974:         R                    = L[1] / (2.0 * PETSC_PI);
975:         ncoords[coordSize++] = PetscSinReal(phi) * r;
976:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
977:         ncoords[coordSize++] = PetscSinReal(theta) * (R + r * PetscCosReal(phi));
978:       } else if (L && (L[0] > 0.0)) {
979:         /* X-periodic */
980:         ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
981:         ncoords[coordSize++] = coords[off + 1];
982:         ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 0]) / L[0]) * (L[0] / (2.0 * PETSC_PI));
983:       } else if (L && (L[1] > 0.0)) {
984:         /* Y-periodic */
985:         ncoords[coordSize++] = coords[off + 0];
986:         ncoords[coordSize++] = PetscSinReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
987:         ncoords[coordSize++] = -PetscCosReal(2.0 * PETSC_PI * PetscRealPart(coords[off + 1]) / L[1]) * (L[1] / (2.0 * PETSC_PI));
988:   #if 0
989:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
990:         PetscReal phi, r, R;
991:         /* Mobius strip */
992:         /* Suppose its an x-z circle, then
993:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
994:            and in that plane we rotate by pi as we go around the circle
995:              \hat r cos(phi/2) + \hat y sin(phi/2) */
996:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
997:         R     = L[0];
998:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
999:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
1000:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
1001:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
1002:   #endif
1003:       } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
1004:     } else {
1005:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off + d];
1006:     }
1007:   }
1008:   if (cutVertexLabel) {
1009:     IS              vertices;
1010:     const PetscInt *verts;
1011:     PetscInt        n;

1013:     DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
1014:     if (vertices) {
1015:       ISGetIndices(vertices, &verts);
1016:       ISGetLocalSize(vertices, &n);
1017:       for (v = 0; v < n; ++v) {
1018:         PetscSectionGetDof(cSection, verts[v], &dof);
1019:         PetscSectionGetOffset(cSection, verts[v], &off);
1020:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off + d] + ((L[d] > 0.) ? L[d] : 0.0);
1021:       }
1022:       ISRestoreIndices(vertices, &verts);
1023:       ISDestroy(&vertices);
1024:     }
1025:   }
1027:   DMLabelDestroy(&cutVertexLabel);
1028:   VecRestoreArray(coordinatesLocal, &coords);
1029:   VecRestoreArray(newcoords, &ncoords);
1030:   PetscObjectSetName((PetscObject)newcoords, "vertices");
1031:   VecScale(newcoords, lengthScale);
1032:   PetscViewerHDF5PushGroup(viewer, "/viz");
1033:   PetscViewerHDF5WriteGroup(viewer, NULL);
1034:   PetscViewerHDF5PopGroup(viewer);
1035:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
1036:   VecView(newcoords, viewer);
1037:   PetscViewerHDF5PopGroup(viewer);
1038:   VecDestroy(&newcoords);
1039:   return 0;
1040: }

1042: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
1043: {
1044:   const char          *topologydm_name;
1045:   const PetscInt      *gpoint;
1046:   PetscInt             numLabels, l;
1047:   DMPlexStorageVersion version;
1048:   char                 group[PETSC_MAX_PATH_LEN];

1050:   PetscViewerHDF5GetDMPlexStorageVersionWriting(viewer, &version);
1051:   ISGetIndices(globalPointNumbers, &gpoint);
1052:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1053:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1054:     PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
1055:   } else {
1056:     PetscStrcpy(group, "/labels");
1057:   }
1058:   PetscViewerHDF5PushGroup(viewer, group);
1059:   DMGetNumLabels(dm, &numLabels);
1060:   for (l = 0; l < numLabels; ++l) {
1061:     DMLabel         label;
1062:     const char     *name;
1063:     IS              valueIS, pvalueIS, globalValueIS;
1064:     const PetscInt *values;
1065:     PetscInt        numValues, v;
1066:     PetscBool       isDepth, output;

1068:     DMGetLabelByNum(dm, l, &label);
1069:     PetscObjectGetName((PetscObject)label, &name);
1070:     DMGetLabelOutput(dm, name, &output);
1071:     PetscStrncmp(name, "depth", 10, &isDepth);
1072:     if (isDepth || !output) continue;
1073:     PetscViewerHDF5PushGroup(viewer, name);
1074:     DMLabelGetValueIS(label, &valueIS);
1075:     /* Must copy to a new IS on the global comm */
1076:     ISGetLocalSize(valueIS, &numValues);
1077:     ISGetIndices(valueIS, &values);
1078:     ISCreateGeneral(PetscObjectComm((PetscObject)dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
1079:     ISRestoreIndices(valueIS, &values);
1080:     ISAllGather(pvalueIS, &globalValueIS);
1081:     ISDestroy(&pvalueIS);
1082:     ISSortRemoveDups(globalValueIS);
1083:     ISGetLocalSize(globalValueIS, &numValues);
1084:     ISGetIndices(globalValueIS, &values);
1085:     for (v = 0; v < numValues; ++v) {
1086:       IS              stratumIS, globalStratumIS;
1087:       const PetscInt *spoints = NULL;
1088:       PetscInt       *gspoints, n = 0, gn, p;
1089:       const char     *iname = "indices";
1090:       char            group[PETSC_MAX_PATH_LEN];

1092:       PetscSNPrintf(group, sizeof(group), "%" PetscInt_FMT, values[v]);
1093:       PetscViewerHDF5PushGroup(viewer, group);
1094:       DMLabelGetStratumIS(label, values[v], &stratumIS);

1096:       if (stratumIS) ISGetLocalSize(stratumIS, &n);
1097:       if (stratumIS) ISGetIndices(stratumIS, &spoints);
1098:       for (gn = 0, p = 0; p < n; ++p)
1099:         if (gpoint[spoints[p]] >= 0) ++gn;
1100:       PetscMalloc1(gn, &gspoints);
1101:       for (gn = 0, p = 0; p < n; ++p)
1102:         if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
1103:       if (stratumIS) ISRestoreIndices(stratumIS, &spoints);
1104:       ISCreateGeneral(PetscObjectComm((PetscObject)dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
1105:       PetscObjectSetName((PetscObject)globalStratumIS, iname);

1107:       ISView(globalStratumIS, viewer);
1108:       ISDestroy(&globalStratumIS);
1109:       ISDestroy(&stratumIS);
1110:       PetscViewerHDF5PopGroup(viewer);
1111:     }
1112:     ISRestoreIndices(globalValueIS, &values);
1113:     ISDestroy(&globalValueIS);
1114:     ISDestroy(&valueIS);
1115:     PetscViewerHDF5PopGroup(viewer);
1116:   }
1117:   ISRestoreIndices(globalPointNumbers, &gpoint);
1118:   PetscViewerHDF5PopGroup(viewer);
1119:   return 0;
1120: }

1122: /* We only write cells and vertices. Does this screw up parallel reading? */
1123: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
1124: {
1125:   IS                globalPointNumbers;
1126:   PetscViewerFormat format;
1127:   PetscBool         viz_geom = PETSC_FALSE, xdmf_topo = PETSC_FALSE, petsc_topo = PETSC_FALSE;

1129:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
1130:   DMPlexCoordinatesView_HDF5_Internal(dm, viewer);

1132:   PetscViewerGetFormat(viewer, &format);
1133:   switch (format) {
1134:   case PETSC_VIEWER_HDF5_VIZ:
1135:     viz_geom  = PETSC_TRUE;
1136:     xdmf_topo = PETSC_TRUE;
1137:     break;
1138:   case PETSC_VIEWER_HDF5_XDMF:
1139:     xdmf_topo = PETSC_TRUE;
1140:     break;
1141:   case PETSC_VIEWER_HDF5_PETSC:
1142:     petsc_topo = PETSC_TRUE;
1143:     break;
1144:   case PETSC_VIEWER_DEFAULT:
1145:   case PETSC_VIEWER_NATIVE:
1146:     viz_geom   = PETSC_TRUE;
1147:     xdmf_topo  = PETSC_TRUE;
1148:     petsc_topo = PETSC_TRUE;
1149:     break;
1150:   default:
1151:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
1152:   }

1154:   if (viz_geom) DMPlexCoordinatesView_HDF5_XDMF_Private(dm, viewer);
1155:   if (xdmf_topo) DMPlexTopologyView_HDF5_XDMF_Private(dm, globalPointNumbers, viewer);
1156:   if (petsc_topo) {
1157:     DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer);
1158:     DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer);
1159:   }

1161:   ISDestroy(&globalPointNumbers);
1162:   return 0;
1163: }

1165: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
1166: {
1167:   MPI_Comm     comm;
1168:   const char  *topologydm_name;
1169:   const char  *sectiondm_name;
1170:   PetscSection gsection;

1172:   PetscObjectGetComm((PetscObject)sectiondm, &comm);
1173:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1174:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1175:   PetscViewerHDF5PushGroup(viewer, "topologies");
1176:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1177:   PetscViewerHDF5PushGroup(viewer, "dms");
1178:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1179:   DMGetGlobalSection(sectiondm, &gsection);
1180:   /* Save raw section */
1181:   PetscSectionView(gsection, viewer);
1182:   /* Save plex wrapper */
1183:   {
1184:     PetscInt        pStart, pEnd, p, n;
1185:     IS              globalPointNumbers;
1186:     const PetscInt *gpoints;
1187:     IS              orderIS;
1188:     PetscInt       *order;

1190:     PetscSectionGetChart(gsection, &pStart, &pEnd);
1191:     DMPlexCreatePointNumbering(dm, &globalPointNumbers);
1192:     ISGetIndices(globalPointNumbers, &gpoints);
1193:     for (p = pStart, n = 0; p < pEnd; ++p)
1194:       if (gpoints[p] >= 0) n++;
1195:     /* "order" is an array of global point numbers.
1196:        When loading, it is used with topology/order array
1197:        to match section points with plex topology points. */
1198:     PetscMalloc1(n, &order);
1199:     for (p = pStart, n = 0; p < pEnd; ++p)
1200:       if (gpoints[p] >= 0) order[n++] = gpoints[p];
1201:     ISRestoreIndices(globalPointNumbers, &gpoints);
1202:     ISDestroy(&globalPointNumbers);
1203:     ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS);
1204:     PetscObjectSetName((PetscObject)orderIS, "order");
1205:     ISView(orderIS, viewer);
1206:     ISDestroy(&orderIS);
1207:   }
1208:   PetscViewerHDF5PopGroup(viewer);
1209:   PetscViewerHDF5PopGroup(viewer);
1210:   PetscViewerHDF5PopGroup(viewer);
1211:   PetscViewerHDF5PopGroup(viewer);
1212:   return 0;
1213: }

1215: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1216: {
1217:   const char *topologydm_name;
1218:   const char *sectiondm_name;
1219:   const char *vec_name;
1220:   PetscInt    bs;

1222:   /* Check consistency */
1223:   {
1224:     PetscSF pointsf, pointsf1;

1226:     DMGetPointSF(dm, &pointsf);
1227:     DMGetPointSF(sectiondm, &pointsf1);
1229:   }
1230:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1231:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1232:   PetscObjectGetName((PetscObject)vec, &vec_name);
1233:   PetscViewerHDF5PushGroup(viewer, "topologies");
1234:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1235:   PetscViewerHDF5PushGroup(viewer, "dms");
1236:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1237:   PetscViewerHDF5PushGroup(viewer, "vecs");
1238:   PetscViewerHDF5PushGroup(viewer, vec_name);
1239:   VecGetBlockSize(vec, &bs);
1240:   PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs);
1241:   VecSetBlockSize(vec, 1);
1242:   /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but,    */
1243:   /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view    */
1244:   /* is set to VecView_Plex, which would save vec in a predefined location.  */
1245:   /* To save vec in where we want, we create a new Vec (temp) with           */
1246:   /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1247:   {
1248:     Vec                temp;
1249:     const PetscScalar *array;
1250:     PetscLayout        map;

1252:     VecCreate(PetscObjectComm((PetscObject)vec), &temp);
1253:     PetscObjectSetName((PetscObject)temp, vec_name);
1254:     VecGetLayout(vec, &map);
1255:     VecSetLayout(temp, map);
1256:     VecSetUp(temp);
1257:     VecGetArrayRead(vec, &array);
1258:     VecPlaceArray(temp, array);
1259:     VecView(temp, viewer);
1260:     VecResetArray(temp);
1261:     VecRestoreArrayRead(vec, &array);
1262:     VecDestroy(&temp);
1263:   }
1264:   VecSetBlockSize(vec, bs);
1265:   PetscViewerHDF5PopGroup(viewer);
1266:   PetscViewerHDF5PopGroup(viewer);
1267:   PetscViewerHDF5PopGroup(viewer);
1268:   PetscViewerHDF5PopGroup(viewer);
1269:   PetscViewerHDF5PopGroup(viewer);
1270:   PetscViewerHDF5PopGroup(viewer);
1271:   return 0;
1272: }

1274: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1275: {
1276:   MPI_Comm     comm;
1277:   const char  *topologydm_name;
1278:   const char  *sectiondm_name;
1279:   const char  *vec_name;
1280:   PetscSection section;
1281:   PetscBool    includesConstraints;
1282:   Vec          gvec;
1283:   PetscInt     m, bs;

1285:   PetscObjectGetComm((PetscObject)dm, &comm);
1286:   /* Check consistency */
1287:   {
1288:     PetscSF pointsf, pointsf1;

1290:     DMGetPointSF(dm, &pointsf);
1291:     DMGetPointSF(sectiondm, &pointsf1);
1293:   }
1294:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1295:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1296:   PetscObjectGetName((PetscObject)vec, &vec_name);
1297:   PetscViewerHDF5PushGroup(viewer, "topologies");
1298:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1299:   PetscViewerHDF5PushGroup(viewer, "dms");
1300:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1301:   PetscViewerHDF5PushGroup(viewer, "vecs");
1302:   PetscViewerHDF5PushGroup(viewer, vec_name);
1303:   VecGetBlockSize(vec, &bs);
1304:   PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *)&bs);
1305:   VecCreate(comm, &gvec);
1306:   PetscObjectSetName((PetscObject)gvec, vec_name);
1307:   DMGetGlobalSection(sectiondm, &section);
1308:   PetscSectionGetIncludesConstraints(section, &includesConstraints);
1309:   if (includesConstraints) PetscSectionGetStorageSize(section, &m);
1310:   else PetscSectionGetConstrainedStorageSize(section, &m);
1311:   VecSetSizes(gvec, m, PETSC_DECIDE);
1312:   VecSetUp(gvec);
1313:   DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec);
1314:   DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec);
1315:   VecView(gvec, viewer);
1316:   VecDestroy(&gvec);
1317:   PetscViewerHDF5PopGroup(viewer);
1318:   PetscViewerHDF5PopGroup(viewer);
1319:   PetscViewerHDF5PopGroup(viewer);
1320:   PetscViewerHDF5PopGroup(viewer);
1321:   PetscViewerHDF5PopGroup(viewer);
1322:   PetscViewerHDF5PopGroup(viewer);
1323:   return 0;
1324: }

1326: struct _n_LoadLabelsCtx {
1327:   MPI_Comm    comm;
1328:   PetscMPIInt rank;
1329:   DM          dm;
1330:   PetscViewer viewer;
1331:   DMLabel     label;
1332:   PetscSF     sfXC;
1333:   PetscLayout layoutX;
1334: };
1335: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;

1337: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1338: {
1339:   PetscNew(ctx);
1340:   PetscObjectReference((PetscObject)((*ctx)->dm = dm));
1341:   PetscObjectReference((PetscObject)((*ctx)->viewer = viewer));
1342:   PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm);
1343:   MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank);
1344:   (*ctx)->sfXC = sfXC;
1345:   if (sfXC) {
1346:     PetscInt nX;

1348:     PetscObjectReference((PetscObject)sfXC);
1349:     PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL);
1350:     PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX);
1351:   }
1352:   return 0;
1353: }

1355: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1356: {
1357:   if (!*ctx) return 0;
1358:   DMDestroy(&(*ctx)->dm);
1359:   PetscViewerDestroy(&(*ctx)->viewer);
1360:   PetscSFDestroy(&(*ctx)->sfXC);
1361:   PetscLayoutDestroy(&(*ctx)->layoutX);
1362:   PetscFree(*ctx);
1363:   return 0;
1364: }

1366: /*
1367:     A: on-disk points
1368:     X: global points [0, NX)
1369:     C: distributed plex points
1370: */
1371: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1372: {
1373:   MPI_Comm        comm    = ctx->comm;
1374:   PetscSF         sfXC    = ctx->sfXC;
1375:   PetscLayout     layoutX = ctx->layoutX;
1376:   PetscSF         sfXA;
1377:   const PetscInt *A_points;
1378:   PetscInt        nX, nC;
1379:   PetscInt        n;

1381:   PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL);
1382:   ISGetLocalSize(stratumIS, &n);
1383:   ISGetIndices(stratumIS, &A_points);
1384:   PetscSFCreate(comm, &sfXA);
1385:   PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points);
1386:   ISCreate(comm, newStratumIS);
1387:   ISSetType(*newStratumIS, ISGENERAL);
1388:   {
1389:     PetscInt   i;
1390:     PetscBool *A_mask, *X_mask, *C_mask;

1392:     PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask);
1393:     for (i = 0; i < n; i++) A_mask[i] = PETSC_TRUE;
1394:     PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1395:     PetscSFReduceEnd(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1396:     PetscSFBcastBegin(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1397:     PetscSFBcastEnd(sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1398:     ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask);
1399:     PetscFree3(A_mask, X_mask, C_mask);
1400:   }
1401:   PetscSFDestroy(&sfXA);
1402:   ISRestoreIndices(stratumIS, &A_points);
1403:   return 0;
1404: }

1406: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1407: {
1408:   LoadLabelsCtx   ctx    = (LoadLabelsCtx)op_data;
1409:   PetscViewer     viewer = ctx->viewer;
1410:   DMLabel         label  = ctx->label;
1411:   MPI_Comm        comm   = ctx->comm;
1412:   IS              stratumIS;
1413:   const PetscInt *ind;
1414:   PetscInt        value, N, i;

1416:   PetscOptionsStringToInt(vname, &value);
1417:   ISCreate(comm, &stratumIS);
1418:   PetscObjectSetName((PetscObject)stratumIS, "indices");
1419:   PetscViewerHDF5PushGroup(viewer, vname); /* labels/<lname>/<vname> */

1421:   if (!ctx->sfXC) {
1422:     /* Force serial load */
1423:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
1424:     PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0);
1425:     PetscLayoutSetSize(stratumIS->map, N);
1426:   }
1427:   ISLoad(stratumIS, viewer);

1429:   if (ctx->sfXC) {
1430:     IS newStratumIS;

1432:     ReadLabelStratumHDF5_Distribute_Private(stratumIS, ctx, &newStratumIS);
1433:     ISDestroy(&stratumIS);
1434:     stratumIS = newStratumIS;
1435:   }

1437:   PetscViewerHDF5PopGroup(viewer);
1438:   ISGetLocalSize(stratumIS, &N);
1439:   ISGetIndices(stratumIS, &ind);
1440:   for (i = 0; i < N; ++i) DMLabelSetValue(label, ind[i], value);
1441:   ISRestoreIndices(stratumIS, &ind);
1442:   ISDestroy(&stratumIS);
1443:   return 0;
1444: }

1446: /* TODO: Fix this code, it is returning PETSc error codes when it should be translating them to herr_t codes */
1447: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1448: {
1449:   LoadLabelsCtx  ctx = (LoadLabelsCtx)op_data;
1450:   DM             dm  = ctx->dm;
1451:   hsize_t        idx = 0;
1453:   PetscBool      flg;
1454:   herr_t         err;

1456:   DMHasLabel(dm, lname, &flg);
1457:   if (flg) DMRemoveLabel(dm, lname, NULL);
1458:   DMCreateLabel(dm, lname);
1459:   if (ierr) return (herr_t)ierr;
1460:   DMGetLabel(dm, lname, &ctx->label);
1461:   if (ierr) return (herr_t)ierr;
1462:   PetscViewerHDF5PushGroup(ctx->viewer, lname);
1463:   if (ierr) return (herr_t)ierr;
1464:   /* Iterate over the label's strata */
1465:   PetscCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1466:   PetscViewerHDF5PopGroup(ctx->viewer);
1467:   if (ierr) return (herr_t)ierr;
1468:   return err;
1469: }

1471: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1472: {
1473:   const char          *topologydm_name;
1474:   LoadLabelsCtx        ctx;
1475:   hsize_t              idx = 0;
1476:   char                 group[PETSC_MAX_PATH_LEN];
1477:   DMPlexStorageVersion version;
1478:   PetscBool            distributed, hasGroup;

1480:   DMPlexIsDistributed(dm, &distributed);
1482:   LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx);
1483:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1484:   PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1485:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1486:     PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
1487:   } else {
1488:     PetscStrcpy(group, "labels");
1489:   }
1490:   PetscViewerHDF5PushGroup(viewer, group);
1491:   PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup);
1492:   if (hasGroup) {
1493:     hid_t fileId, groupId;

1495:     PetscViewerHDF5OpenGroup(viewer, NULL, &fileId, &groupId);
1496:     /* Iterate over labels */
1497:     PetscCallHDF5(H5Literate, (groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1498:     PetscCallHDF5(H5Gclose, (groupId));
1499:   }
1500:   PetscViewerHDF5PopGroup(viewer);
1501:   LoadLabelsCtxDestroy(&ctx);
1502:   return 0;
1503: }

1505: static PetscErrorCode DMPlexDistributionLoad_HDF5_Private(DM dm, PetscViewer viewer, PetscSF sf, PetscSF *distsf, DM *distdm)
1506: {
1507:   MPI_Comm        comm;
1508:   PetscMPIInt     size, rank, dist_size;
1509:   const char     *distribution_name;
1510:   PetscInt        p, lsize;
1511:   IS              chartSizesIS, ownersIS, gpointsIS;
1512:   const PetscInt *chartSize, *owners, *gpoints;
1513:   PetscLayout     layout;
1514:   PetscBool       has;

1516:   *distsf = NULL;
1517:   *distdm = NULL;
1518:   PetscObjectGetComm((PetscObject)dm, &comm);
1519:   MPI_Comm_size(comm, &size);
1520:   MPI_Comm_rank(comm, &rank);
1521:   DMPlexDistributionGetName(dm, &distribution_name);
1522:   if (!distribution_name) return 0;
1523:   PetscViewerHDF5HasGroup(viewer, NULL, &has);
1524:   if (!has) {
1525:     char *full_group;

1527:     PetscViewerHDF5GetGroup(viewer, NULL, &full_group);
1529:   }
1530:   PetscViewerHDF5ReadAttribute(viewer, NULL, "comm_size", PETSC_INT, NULL, (void *)&dist_size);
1532:   ISCreate(comm, &chartSizesIS);
1533:   PetscObjectSetName((PetscObject)chartSizesIS, "chart_sizes");
1534:   ISCreate(comm, &ownersIS);
1535:   PetscObjectSetName((PetscObject)ownersIS, "owners");
1536:   ISCreate(comm, &gpointsIS);
1537:   PetscObjectSetName((PetscObject)gpointsIS, "global_point_numbers");
1538:   PetscLayoutSetLocalSize(chartSizesIS->map, 1);
1539:   ISLoad(chartSizesIS, viewer);
1540:   ISGetIndices(chartSizesIS, &chartSize);
1541:   PetscLayoutSetLocalSize(ownersIS->map, *chartSize);
1542:   PetscLayoutSetLocalSize(gpointsIS->map, *chartSize);
1543:   ISLoad(ownersIS, viewer);
1544:   ISLoad(gpointsIS, viewer);
1545:   ISGetIndices(ownersIS, &owners);
1546:   ISGetIndices(gpointsIS, &gpoints);
1547:   PetscSFCreate(comm, distsf);
1548:   PetscSFSetFromOptions(*distsf);
1549:   PetscLayoutCreate(comm, &layout);
1550:   PetscSFGetGraph(sf, &lsize, NULL, NULL, NULL);
1551:   PetscLayoutSetLocalSize(layout, lsize);
1552:   PetscLayoutSetBlockSize(layout, 1);
1553:   PetscLayoutSetUp(layout);
1554:   PetscSFSetGraphLayout(*distsf, layout, *chartSize, NULL, PETSC_OWN_POINTER, gpoints);
1555:   PetscLayoutDestroy(&layout);
1556:   /* Migrate DM */
1557:   {
1558:     PetscInt     pStart, pEnd;
1559:     PetscSFNode *buffer0, *buffer1, *buffer2;

1561:     DMPlexGetChart(dm, &pStart, &pEnd);
1562:     PetscMalloc2(pEnd - pStart, &buffer0, lsize, &buffer1);
1563:     PetscMalloc1(*chartSize, &buffer2);
1564:     {
1565:       PetscSF            workPointSF;
1566:       PetscInt           workNroots, workNleaves;
1567:       const PetscInt    *workIlocal;
1568:       const PetscSFNode *workIremote;

1570:       for (p = pStart; p < pEnd; ++p) {
1571:         buffer0[p - pStart].rank  = rank;
1572:         buffer0[p - pStart].index = p - pStart;
1573:       }
1574:       DMGetPointSF(dm, &workPointSF);
1575:       PetscSFGetGraph(workPointSF, &workNroots, &workNleaves, &workIlocal, &workIremote);
1576:       for (p = 0; p < workNleaves; ++p) {
1577:         PetscInt workIlocalp = (workIlocal ? workIlocal[p] : p);

1579:         buffer0[workIlocalp].rank = -1;
1580:       }
1581:     }
1582:     for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1583:     for (p = 0; p < *chartSize; ++p) buffer2[p].rank = -1;
1584:     PetscSFReduceBegin(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1585:     PetscSFReduceEnd(sf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1586:     PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE);
1587:     PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer2, MPI_REPLACE);
1588:     if (PetscDefined(USE_DEBUG)) {
1589:       for (p = 0; p < *chartSize; ++p) {
1591:       }
1592:     }
1593:     PetscFree2(buffer0, buffer1);
1594:     DMCreate(comm, distdm);
1595:     DMSetType(*distdm, DMPLEX);
1596:     PetscObjectSetName((PetscObject)(*distdm), ((PetscObject)dm)->name);
1597:     DMPlexDistributionSetName(*distdm, distribution_name);
1598:     {
1599:       PetscSF migrationSF;

1601:       PetscSFCreate(comm, &migrationSF);
1602:       PetscSFSetFromOptions(migrationSF);
1603:       PetscSFSetGraph(migrationSF, pEnd - pStart, *chartSize, NULL, PETSC_OWN_POINTER, buffer2, PETSC_OWN_POINTER);
1604:       PetscSFSetUp(migrationSF);
1605:       DMPlexMigrate(dm, migrationSF, *distdm);
1606:       PetscSFDestroy(&migrationSF);
1607:     }
1608:   }
1609:   /* Set pointSF */
1610:   {
1611:     PetscSF      pointSF;
1612:     PetscInt    *ilocal, nleaves, q;
1613:     PetscSFNode *iremote, *buffer0, *buffer1;

1615:     PetscMalloc2(*chartSize, &buffer0, lsize, &buffer1);
1616:     for (p = 0, nleaves = 0; p < *chartSize; ++p) {
1617:       if (owners[p] == rank) {
1618:         buffer0[p].rank = rank;
1619:       } else {
1620:         buffer0[p].rank = -1;
1621:         nleaves++;
1622:       }
1623:       buffer0[p].index = p;
1624:     }
1625:     for (p = 0; p < lsize; ++p) buffer1[p].rank = -1;
1626:     PetscSFReduceBegin(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1627:     PetscSFReduceEnd(*distsf, MPIU_2INT, buffer0, buffer1, MPI_MAXLOC);
1628:     for (p = 0; p < *chartSize; ++p) buffer0[p].rank = -1;
1629:     PetscSFBcastBegin(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE);
1630:     PetscSFBcastEnd(*distsf, MPIU_2INT, buffer1, buffer0, MPI_REPLACE);
1631:     if (PetscDefined(USE_DEBUG)) {
1633:     }
1634:     PetscMalloc1(nleaves, &ilocal);
1635:     PetscMalloc1(nleaves, &iremote);
1636:     for (p = 0, q = 0; p < *chartSize; ++p) {
1637:       if (buffer0[p].rank != rank) {
1638:         ilocal[q]        = p;
1639:         iremote[q].rank  = buffer0[p].rank;
1640:         iremote[q].index = buffer0[p].index;
1641:         q++;
1642:       }
1643:     }
1645:     PetscFree2(buffer0, buffer1);
1646:     PetscSFCreate(comm, &pointSF);
1647:     PetscSFSetFromOptions(pointSF);
1648:     PetscSFSetGraph(pointSF, *chartSize, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1649:     DMSetPointSF(*distdm, pointSF);
1650:     {
1651:       DM cdm;

1653:       DMGetCoordinateDM(*distdm, &cdm);
1654:       DMSetPointSF(cdm, pointSF);
1655:     }
1656:     PetscSFDestroy(&pointSF);
1657:   }
1658:   ISRestoreIndices(chartSizesIS, &chartSize);
1659:   ISRestoreIndices(ownersIS, &owners);
1660:   ISRestoreIndices(gpointsIS, &gpoints);
1661:   ISDestroy(&chartSizesIS);
1662:   ISDestroy(&ownersIS);
1663:   ISDestroy(&gpointsIS);
1664:   /* Record that overlap has been manually created.               */
1665:   /* This is to pass `DMPlexCheckPointSF()`, which checks that    */
1666:   /* pointSF does not contain cells in the leaves if overlap = 0. */
1667:   DMPlexSetOverlap_Plex(*distdm, NULL, DMPLEX_OVERLAP_MANUAL);
1668:   DMPlexDistributeSetDefault(*distdm, PETSC_FALSE);
1669:   DMPlexReorderSetDefault(*distdm, DMPLEX_REORDER_DEFAULT_FALSE);
1670:   return 0;
1671: }

1673: static PetscErrorCode DMPlexTopologyLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer, PetscSF *sf)
1674: {
1675:   MPI_Comm        comm;
1676:   const char     *pointsName, *coneSizesName, *conesName, *orientationsName;
1677:   IS              pointsIS, coneSizesIS, conesIS, orientationsIS;
1678:   const PetscInt *points, *coneSizes, *cones, *orientations;
1679:   PetscInt       *cone, *ornt;
1680:   PetscInt        dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1681:   PetscMPIInt     size, rank;

1683:   pointsName       = "order";
1684:   coneSizesName    = "cones";
1685:   conesName        = "cells";
1686:   orientationsName = "orientation";
1687:   PetscObjectGetComm((PetscObject)dm, &comm);
1688:   MPI_Comm_size(comm, &size);
1689:   MPI_Comm_rank(comm, &rank);
1690:   ISCreate(comm, &pointsIS);
1691:   PetscObjectSetName((PetscObject)pointsIS, pointsName);
1692:   ISCreate(comm, &coneSizesIS);
1693:   PetscObjectSetName((PetscObject)coneSizesIS, coneSizesName);
1694:   ISCreate(comm, &conesIS);
1695:   PetscObjectSetName((PetscObject)conesIS, conesName);
1696:   ISCreate(comm, &orientationsIS);
1697:   PetscObjectSetName((PetscObject)orientationsIS, orientationsName);
1698:   PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject)conesIS, "cell_dim", PETSC_INT, NULL, &dim);
1699:   DMSetDimension(dm, dim);
1700:   {
1701:     /* Force serial load */
1702:     PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np);
1703:     PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0);
1704:     PetscLayoutSetSize(pointsIS->map, Np);
1705:     pEnd = rank == 0 ? Np : 0;
1706:     PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np);
1707:     PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0);
1708:     PetscLayoutSetSize(coneSizesIS->map, Np);
1709:     PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N);
1710:     PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0);
1711:     PetscLayoutSetSize(conesIS->map, N);
1712:     PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N);
1713:     PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0);
1714:     PetscLayoutSetSize(orientationsIS->map, N);
1715:   }
1716:   ISLoad(pointsIS, viewer);
1717:   ISLoad(coneSizesIS, viewer);
1718:   ISLoad(conesIS, viewer);
1719:   ISLoad(orientationsIS, viewer);
1720:   /* Create Plex */
1721:   DMPlexSetChart(dm, 0, pEnd);
1722:   ISGetIndices(pointsIS, &points);
1723:   ISGetIndices(coneSizesIS, &coneSizes);
1724:   for (p = 0; p < pEnd; ++p) {
1725:     DMPlexSetConeSize(dm, points[p], coneSizes[p]);
1726:     maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
1727:   }
1728:   DMSetUp(dm);
1729:   ISGetIndices(conesIS, &cones);
1730:   ISGetIndices(orientationsIS, &orientations);
1731:   PetscMalloc2(maxConeSize, &cone, maxConeSize, &ornt);
1732:   for (p = 0, q = 0; p < pEnd; ++p) {
1733:     for (c = 0; c < coneSizes[p]; ++c, ++q) {
1734:       cone[c] = cones[q];
1735:       ornt[c] = orientations[q];
1736:     }
1737:     DMPlexSetCone(dm, points[p], cone);
1738:     DMPlexSetConeOrientation(dm, points[p], ornt);
1739:   }
1740:   PetscFree2(cone, ornt);
1741:   /* Create global section migration SF */
1742:   if (sf) {
1743:     PetscLayout layout;
1744:     PetscInt   *globalIndices;

1746:     PetscMalloc1(pEnd, &globalIndices);
1747:     /* plex point == globalPointNumber in this case */
1748:     for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1749:     PetscLayoutCreate(comm, &layout);
1750:     PetscLayoutSetSize(layout, Np);
1751:     PetscLayoutSetBlockSize(layout, 1);
1752:     PetscLayoutSetUp(layout);
1753:     PetscSFCreate(comm, sf);
1754:     PetscSFSetFromOptions(*sf);
1755:     PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices);
1756:     PetscLayoutDestroy(&layout);
1757:     PetscFree(globalIndices);
1758:   }
1759:   /* Clean-up */
1760:   ISRestoreIndices(pointsIS, &points);
1761:   ISRestoreIndices(coneSizesIS, &coneSizes);
1762:   ISRestoreIndices(conesIS, &cones);
1763:   ISRestoreIndices(orientationsIS, &orientations);
1764:   ISDestroy(&pointsIS);
1765:   ISDestroy(&coneSizesIS);
1766:   ISDestroy(&conesIS);
1767:   ISDestroy(&orientationsIS);
1768:   /* Fill in the rest of the topology structure */
1769:   DMPlexSymmetrize(dm);
1770:   DMPlexStratify(dm);
1771:   return 0;
1772: }

1774: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sfXC)
1775: {
1776:   DMPlexStorageVersion version;
1777:   const char          *topologydm_name;
1778:   char                 group[PETSC_MAX_PATH_LEN];
1779:   PetscSF              sfwork = NULL;

1781:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1782:   PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1783:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1784:     PetscSNPrintf(group, sizeof(group), "topologies/%s", topologydm_name);
1785:   } else {
1786:     PetscStrcpy(group, "/");
1787:   }
1788:   PetscViewerHDF5PushGroup(viewer, group);

1790:   PetscViewerHDF5PushGroup(viewer, "topology");
1791:   DMPlexTopologyLoad_HDF5_Legacy_Private(dm, viewer, &sfwork);
1792:   PetscViewerHDF5PopGroup(viewer); /* "topology" */

1794:   if (DMPlexStorageVersionGE(version, 2, 0, 0)) {
1795:     DM          distdm;
1796:     PetscSF     distsf;
1797:     const char *distribution_name;
1798:     PetscBool   exists;

1800:     DMPlexDistributionGetName(dm, &distribution_name);
1801:     /* this group is guaranteed to be present since DMPlexStorageVersion 2.1.0 */
1802:     PetscViewerHDF5PushGroup(viewer, "distributions");
1803:     PetscViewerHDF5HasGroup(viewer, NULL, &exists);
1804:     if (exists) {
1805:       PetscViewerHDF5PushGroup(viewer, distribution_name);
1806:       DMPlexDistributionLoad_HDF5_Private(dm, viewer, sfwork, &distsf, &distdm);
1807:       if (distdm) {
1808:         DMPlexReplace_Internal(dm, &distdm);
1809:         PetscSFDestroy(&sfwork);
1810:         sfwork = distsf;
1811:       }
1812:       PetscViewerHDF5PopGroup(viewer); /* distribution_name */
1813:     }
1814:     PetscViewerHDF5PopGroup(viewer); /* "distributions" */
1815:   }
1816:   if (sfXC) {
1817:     *sfXC = sfwork;
1818:   } else {
1819:     PetscSFDestroy(&sfwork);
1820:   }

1822:   PetscViewerHDF5PopGroup(viewer);
1823:   return 0;
1824: }

1826: /* If the file is old, it not only has different path to the coordinates, but   */
1827: /* does not contain coordinateDMs, so must fall back to the old implementation. */
1828: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1829: {
1830:   PetscSection coordSection;
1831:   Vec          coordinates;
1832:   PetscReal    lengthScale;
1833:   PetscInt     spatialDim, N, numVertices, vStart, vEnd, v;
1834:   PetscMPIInt  rank;

1836:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1837:   /* Read geometry */
1838:   PetscViewerHDF5PushGroup(viewer, "/geometry");
1839:   VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);
1840:   PetscObjectSetName((PetscObject)coordinates, "vertices");
1841:   {
1842:     /* Force serial load */
1843:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
1844:     VecSetSizes(coordinates, rank == 0 ? N : 0, N);
1845:     VecSetBlockSize(coordinates, spatialDim);
1846:   }
1847:   VecLoad(coordinates, viewer);
1848:   PetscViewerHDF5PopGroup(viewer);
1849:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1850:   VecScale(coordinates, 1.0 / lengthScale);
1851:   VecGetLocalSize(coordinates, &numVertices);
1852:   VecGetBlockSize(coordinates, &spatialDim);
1853:   numVertices /= spatialDim;
1854:   /* Create coordinates */
1855:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1857:   DMGetCoordinateSection(dm, &coordSection);
1858:   PetscSectionSetNumFields(coordSection, 1);
1859:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
1860:   PetscSectionSetChart(coordSection, vStart, vEnd);
1861:   for (v = vStart; v < vEnd; ++v) {
1862:     PetscSectionSetDof(coordSection, v, spatialDim);
1863:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1864:   }
1865:   PetscSectionSetUp(coordSection);
1866:   DMSetCoordinates(dm, coordinates);
1867:   VecDestroy(&coordinates);
1868:   return 0;
1869: }

1871: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1872: {
1873:   DMPlexStorageVersion version;
1874:   DM                   cdm;
1875:   Vec                  coords;
1876:   PetscInt             blockSize;
1877:   PetscReal            lengthScale;
1878:   PetscSF              lsf;
1879:   const char          *topologydm_name;
1880:   char                *coordinatedm_name, *coordinates_name;

1882:   PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1883:   if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
1884:     DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1885:     return 0;
1886:   }
1887:   /* else: since DMPlexStorageVersion 2.0.0 */
1889:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1890:   PetscViewerHDF5PushGroup(viewer, "topologies");
1891:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1892:   PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, NULL, &coordinatedm_name);
1893:   PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, NULL, &coordinates_name);
1894:   PetscViewerHDF5PopGroup(viewer);
1895:   PetscViewerHDF5PopGroup(viewer);
1896:   DMGetCoordinateDM(dm, &cdm);
1897:   PetscObjectSetName((PetscObject)cdm, coordinatedm_name);
1898:   PetscFree(coordinatedm_name);
1899:   /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
1900:   DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf);
1901:   DMCreateLocalVector(cdm, &coords);
1902:   PetscObjectSetName((PetscObject)coords, coordinates_name);
1903:   PetscFree(coordinates_name);
1904:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
1905:   DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords);
1906:   PetscViewerPopFormat(viewer);
1907:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1908:   VecScale(coords, 1.0 / lengthScale);
1909:   DMSetCoordinatesLocal(dm, coords);
1910:   VecGetBlockSize(coords, &blockSize);
1911:   DMSetCoordinateDim(dm, blockSize);
1912:   VecDestroy(&coords);
1913:   PetscSFDestroy(&lsf);
1914:   return 0;
1915: }

1917: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1918: {
1919:   DMPlexStorageVersion version;

1921:   PetscViewerHDF5GetDMPlexStorageVersionReading(viewer, &version);
1922:   if (!DMPlexStorageVersionGE(version, 2, 0, 0)) {
1923:     DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL);
1924:     DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL);
1925:     DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1926:   } else {
1927:     PetscSF sfXC;

1929:     /* since DMPlexStorageVersion 2.0.0 */
1930:     DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC);
1931:     DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC);
1932:     DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC);
1933:     PetscSFDestroy(&sfXC);
1934:   }
1935:   return 0;
1936: }

1938: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1939: {
1940:   MPI_Comm  comm;
1941:   PetscInt  pStart, pEnd, p, m;
1942:   PetscInt *goffs, *ilocal;
1943:   PetscBool rootIncludeConstraints, leafIncludeConstraints;

1945:   PetscObjectGetComm((PetscObject)leafSection, &comm);
1946:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1947:   PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints);
1948:   PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints);
1949:   if (rootIncludeConstraints && leafIncludeConstraints) PetscSectionGetStorageSize(leafSection, &m);
1950:   else PetscSectionGetConstrainedStorageSize(leafSection, &m);
1951:   PetscMalloc1(m, &ilocal);
1952:   PetscMalloc1(m, &goffs);
1953:   /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
1954:   /* for the top-level section (not for each field), so one must have   */
1955:   /* rootSection->pointMajor == PETSC_TRUE.                             */
1957:   /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
1959:   for (p = pStart, m = 0; p < pEnd; ++p) {
1960:     PetscInt        dof, cdof, i, j, off, goff;
1961:     const PetscInt *cinds;

1963:     PetscSectionGetDof(leafSection, p, &dof);
1964:     if (dof < 0) continue;
1965:     goff = globalOffsets[p - pStart];
1966:     PetscSectionGetOffset(leafSection, p, &off);
1967:     PetscSectionGetConstraintDof(leafSection, p, &cdof);
1968:     PetscSectionGetConstraintIndices(leafSection, p, &cinds);
1969:     for (i = 0, j = 0; i < dof; ++i) {
1970:       PetscBool constrained = (PetscBool)(j < cdof && i == cinds[j]);

1972:       if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {
1973:         ilocal[m]  = off++;
1974:         goffs[m++] = goff++;
1975:       } else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
1976:       else if (!leafIncludeConstraints && rootIncludeConstraints) ++goff;
1977:       if (constrained) ++j;
1978:     }
1979:   }
1980:   PetscSFCreate(comm, sectionSF);
1981:   PetscSFSetFromOptions(*sectionSF);
1982:   PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs);
1983:   PetscFree(goffs);
1984:   return 0;
1985: }

1987: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
1988: {
1989:   MPI_Comm     comm;
1990:   PetscMPIInt  size, rank;
1991:   const char  *topologydm_name;
1992:   const char  *sectiondm_name;
1993:   PetscSection sectionA, sectionB;
1994:   PetscInt     nX, n, i;
1995:   PetscSF      sfAB;

1997:   PetscObjectGetComm((PetscObject)dm, &comm);
1998:   MPI_Comm_size(comm, &size);
1999:   MPI_Comm_rank(comm, &rank);
2000:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
2001:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
2002:   PetscViewerHDF5PushGroup(viewer, "topologies");
2003:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
2004:   PetscViewerHDF5PushGroup(viewer, "dms");
2005:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
2006:   /* A: on-disk points                        */
2007:   /* X: list of global point numbers, [0, NX) */
2008:   /* B: plex points                           */
2009:   /* Load raw section (sectionA)              */
2010:   PetscSectionCreate(comm, &sectionA);
2011:   PetscSectionLoad(sectionA, viewer);
2012:   PetscSectionGetChart(sectionA, NULL, &n);
2013:   /* Create sfAB: A -> B */
2014:   #if defined(PETSC_USE_DEBUG)
2015:   {
2016:     PetscInt N, N1;

2018:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1);
2019:     MPI_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm);
2021:   }
2022:   #endif
2023:   {
2024:     IS              orderIS;
2025:     const PetscInt *gpoints;
2026:     PetscSF         sfXA, sfAX;
2027:     PetscLayout     layout;
2028:     PetscSFNode    *owners, *buffer;
2029:     PetscInt        nleaves;
2030:     PetscInt       *ilocal;
2031:     PetscSFNode    *iremote;

2033:     /* Create sfAX: A -> X */
2034:     ISCreate(comm, &orderIS);
2035:     PetscObjectSetName((PetscObject)orderIS, "order");
2036:     PetscLayoutSetLocalSize(orderIS->map, n);
2037:     ISLoad(orderIS, viewer);
2038:     PetscLayoutCreate(comm, &layout);
2039:     PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL);
2040:     PetscLayoutSetLocalSize(layout, nX);
2041:     PetscLayoutSetBlockSize(layout, 1);
2042:     PetscLayoutSetUp(layout);
2043:     PetscSFCreate(comm, &sfXA);
2044:     ISGetIndices(orderIS, &gpoints);
2045:     PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints);
2046:     ISRestoreIndices(orderIS, &gpoints);
2047:     ISDestroy(&orderIS);
2048:     PetscLayoutDestroy(&layout);
2049:     PetscMalloc1(n, &owners);
2050:     PetscMalloc1(nX, &buffer);
2051:     for (i = 0; i < n; ++i) {
2052:       owners[i].rank  = rank;
2053:       owners[i].index = i;
2054:     }
2055:     for (i = 0; i < nX; ++i) {
2056:       buffer[i].rank  = -1;
2057:       buffer[i].index = -1;
2058:     }
2059:     PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
2060:     PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
2061:     PetscSFDestroy(&sfXA);
2062:     PetscFree(owners);
2063:     for (i = 0, nleaves = 0; i < nX; ++i)
2064:       if (buffer[i].rank >= 0) nleaves++;
2065:     PetscMalloc1(nleaves, &ilocal);
2066:     PetscMalloc1(nleaves, &iremote);
2067:     for (i = 0, nleaves = 0; i < nX; ++i) {
2068:       if (buffer[i].rank >= 0) {
2069:         ilocal[nleaves]        = i;
2070:         iremote[nleaves].rank  = buffer[i].rank;
2071:         iremote[nleaves].index = buffer[i].index;
2072:         nleaves++;
2073:       }
2074:     }
2075:     PetscFree(buffer);
2076:     PetscSFCreate(comm, &sfAX);
2077:     PetscSFSetFromOptions(sfAX);
2078:     PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
2079:     PetscSFCompose(sfAX, sfXB, &sfAB);
2080:     PetscSFDestroy(&sfAX);
2081:   }
2082:   PetscViewerHDF5PopGroup(viewer);
2083:   PetscViewerHDF5PopGroup(viewer);
2084:   PetscViewerHDF5PopGroup(viewer);
2085:   PetscViewerHDF5PopGroup(viewer);
2086:   /* Create plex section (sectionB) */
2087:   DMGetLocalSection(sectiondm, &sectionB);
2088:   if (lsf || gsf) {
2089:     PetscLayout layout;
2090:     PetscInt    M, m;
2091:     PetscInt   *offsetsA;
2092:     PetscBool   includesConstraintsA;

2094:     PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB);
2095:     PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA);
2096:     if (includesConstraintsA) PetscSectionGetStorageSize(sectionA, &m);
2097:     else PetscSectionGetConstrainedStorageSize(sectionA, &m);
2098:     MPI_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm);
2099:     PetscLayoutCreate(comm, &layout);
2100:     PetscLayoutSetSize(layout, M);
2101:     PetscLayoutSetUp(layout);
2102:     if (lsf) {
2103:       PetscSF lsfABdata;

2105:       DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata);
2106:       *lsf = lsfABdata;
2107:     }
2108:     if (gsf) {
2109:       PetscSection gsectionB, gsectionB1;
2110:       PetscBool    includesConstraintsB;
2111:       PetscSF      gsfABdata, pointsf;

2113:       DMGetGlobalSection(sectiondm, &gsectionB1);
2114:       PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB);
2115:       DMGetPointSF(sectiondm, &pointsf);
2116:       PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB);
2117:       DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata);
2118:       PetscSectionDestroy(&gsectionB);
2119:       *gsf = gsfABdata;
2120:     }
2121:     PetscLayoutDestroy(&layout);
2122:     PetscFree(offsetsA);
2123:   } else {
2124:     PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB);
2125:   }
2126:   PetscSFDestroy(&sfAB);
2127:   PetscSectionDestroy(&sectionA);
2128:   return 0;
2129: }

2131: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
2132: {
2133:   MPI_Comm           comm;
2134:   const char        *topologydm_name;
2135:   const char        *sectiondm_name;
2136:   const char        *vec_name;
2137:   Vec                vecA;
2138:   PetscInt           mA, m, bs;
2139:   const PetscInt    *ilocal;
2140:   const PetscScalar *src;
2141:   PetscScalar       *dest;

2143:   PetscObjectGetComm((PetscObject)dm, &comm);
2144:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
2145:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
2146:   PetscObjectGetName((PetscObject)vec, &vec_name);
2147:   PetscViewerHDF5PushGroup(viewer, "topologies");
2148:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
2149:   PetscViewerHDF5PushGroup(viewer, "dms");
2150:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
2151:   PetscViewerHDF5PushGroup(viewer, "vecs");
2152:   PetscViewerHDF5PushGroup(viewer, vec_name);
2153:   VecCreate(comm, &vecA);
2154:   PetscObjectSetName((PetscObject)vecA, vec_name);
2155:   PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL);
2156:   /* Check consistency */
2157:   {
2158:     PetscSF  pointsf, pointsf1;
2159:     PetscInt m1, i, j;

2161:     DMGetPointSF(dm, &pointsf);
2162:     DMGetPointSF(sectiondm, &pointsf1);
2164:   #if defined(PETSC_USE_DEBUG)
2165:     {
2166:       PetscInt MA, MA1;

2168:       MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm);
2169:       PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1);
2171:     }
2172:   #endif
2173:     VecGetLocalSize(vec, &m1);
2175:     for (i = 0; i < m; ++i) {
2176:       j = ilocal ? ilocal[i] : i;
2178:     }
2179:   }
2180:   VecSetSizes(vecA, mA, PETSC_DECIDE);
2181:   VecLoad(vecA, viewer);
2182:   VecGetArrayRead(vecA, &src);
2183:   VecGetArray(vec, &dest);
2184:   PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
2185:   PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
2186:   VecRestoreArray(vec, &dest);
2187:   VecRestoreArrayRead(vecA, &src);
2188:   VecDestroy(&vecA);
2189:   PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *)&bs);
2190:   VecSetBlockSize(vec, bs);
2191:   PetscViewerHDF5PopGroup(viewer);
2192:   PetscViewerHDF5PopGroup(viewer);
2193:   PetscViewerHDF5PopGroup(viewer);
2194:   PetscViewerHDF5PopGroup(viewer);
2195:   PetscViewerHDF5PopGroup(viewer);
2196:   PetscViewerHDF5PopGroup(viewer);
2197:   return 0;
2198: }
2199: #endif