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)

  9: typedef struct DMPlexStorageVersion {
 10:   PetscInt major, minor, subminor;
 11: } DMPlexStorageVersion;

 13: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

 15: static PetscErrorCode DMPlexStorageVersionParseString_Private(DM dm, const char str[], DMPlexStorageVersion *v)
 16: {
 17:   PetscToken      t;
 18:   char           *ts;
 19:   PetscInt        i;
 20:   PetscInt        ti[3];

 22:   PetscTokenCreate(str, '.', &t);
 23:   for (i=0; i<3; i++) {
 24:     PetscTokenFind(t, &ts);
 26:     PetscOptionsStringToInt(ts, &ti[i]);
 27:   }
 28:   PetscTokenFind(t, &ts);
 30:   PetscTokenDestroy(&t);
 31:   v->major    = ti[0];
 32:   v->minor    = ti[1];
 33:   v->subminor = ti[2];
 34:   return 0;
 35: }

 37: static PetscErrorCode DMPlexStorageVersionSetUpWriting_Private(DM dm, PetscViewer viewer, DMPlexStorageVersion *version)
 38: {
 39:   const char      ATTR_NAME[] = "dmplex_storage_version";
 40:   PetscBool       fileHasVersion;
 41:   char            optVersion[16], fileVersion[16];
 42:   PetscErrorCode  ierr;

 44:   PetscStrcpy(fileVersion, DMPLEX_STORAGE_VERSION_STABLE);
 45:   PetscViewerHDF5HasAttribute(viewer, NULL, ATTR_NAME, &fileHasVersion);
 46:   if (fileHasVersion) {
 47:     char *tmp;

 49:     PetscViewerHDF5ReadAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, NULL, &tmp);
 50:     PetscStrcpy(fileVersion, tmp);
 51:     PetscFree(tmp);
 52:   }
 53:   PetscStrcpy(optVersion, fileVersion);
 54:   PetscOptionsBegin(PetscObjectComm((PetscObject)dm),((PetscObject)dm)->prefix,"DMPlex HDF5 Viewer Options","PetscViewer");
 55:   PetscOptionsString("-dm_plex_view_hdf5_storage_version","DMPlex HDF5 viewer storage version",NULL,optVersion,optVersion,sizeof(optVersion),NULL);
 56:   PetscOptionsEnd();
 57:   if (!fileHasVersion) {
 58:     PetscViewerHDF5WriteAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, optVersion);
 59:   } else {
 60:     PetscBool flg;

 62:     PetscStrcmp(fileVersion, optVersion, &flg);
 64:   }
 65:   PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);
 66:   DMPlexStorageVersionParseString_Private(dm, optVersion, version);
 67:   return 0;
 68: }

 70: static PetscErrorCode DMPlexStorageVersionGet_Private(DM dm, PetscViewer viewer, DMPlexStorageVersion *version)
 71: {
 72:   const char      ATTR_NAME[]       = "dmplex_storage_version";
 73:   char           *defaultVersion;
 74:   char           *versionString;

 76:   //TODO string HDF5 attribute handling is terrible and should be redesigned
 77:   PetscStrallocpy("1.0.0", &defaultVersion);
 78:   PetscViewerHDF5ReadAttribute(viewer, NULL, ATTR_NAME, PETSC_STRING, &defaultVersion, &versionString);
 79:   DMPlexStorageVersionParseString_Private(dm, versionString, version);
 80:   PetscFree(versionString);
 81:   PetscFree(defaultVersion);
 82:   return 0;
 83: }

 85: static PetscErrorCode DMPlexGetHDF5Name_Private(DM dm, const char *name[])
 86: {
 87:   if (((PetscObject)dm)->name) {
 88:     PetscObjectGetName((PetscObject)dm, name);
 89:   } else {
 90:     *name = "plex";
 91:   }
 92:   return 0;
 93: }

 95: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
 96: {
 97:   Vec            stamp;
 98:   PetscMPIInt    rank;

100:   if (seqnum < 0) return 0;
101:   MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
102:   VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
103:   VecSetBlockSize(stamp, 1);
104:   PetscObjectSetName((PetscObject) stamp, seqname);
105:   if (rank == 0) {
106:     PetscReal timeScale;
107:     PetscBool istime;

109:     PetscStrncmp(seqname, "time", 5, &istime);
110:     if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); value *= timeScale;}
111:     VecSetValue(stamp, 0, value, INSERT_VALUES);
112:   }
113:   VecAssemblyBegin(stamp);
114:   VecAssemblyEnd(stamp);
115:   PetscViewerHDF5PushGroup(viewer, "/");
116:   PetscViewerHDF5PushTimestepping(viewer);
117:   PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
118:   VecView(stamp, viewer);
119:   PetscViewerHDF5PopTimestepping(viewer);
120:   PetscViewerHDF5PopGroup(viewer);
121:   VecDestroy(&stamp);
122:   return 0;
123: }

125: PetscErrorCode DMSequenceLoad_HDF5_Internal(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
126: {
127:   Vec            stamp;
128:   PetscMPIInt    rank;

130:   if (seqnum < 0) return 0;
131:   MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
132:   VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
133:   VecSetBlockSize(stamp, 1);
134:   PetscObjectSetName((PetscObject) stamp, seqname);
135:   PetscViewerHDF5PushGroup(viewer, "/");
136:   PetscViewerHDF5PushTimestepping(viewer);
137:   PetscViewerHDF5SetTimestep(viewer, seqnum);  /* seqnum < 0 jumps out above */
138:   VecLoad(stamp, viewer);
139:   PetscViewerHDF5PopTimestepping(viewer);
140:   PetscViewerHDF5PopGroup(viewer);
141:   if (rank == 0) {
142:     const PetscScalar *a;
143:     PetscReal timeScale;
144:     PetscBool istime;

146:     VecGetArrayRead(stamp, &a);
147:     *value = a[0];
148:     VecRestoreArrayRead(stamp, &a);
149:     PetscStrncmp(seqname, "time", 5, &istime);
150:     if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
151:   }
152:   VecDestroy(&stamp);
153:   return 0;
154: }

156: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
157: {
158:   IS              cutcells = NULL;
159:   const PetscInt *cutc;
160:   PetscInt        cellHeight, vStart, vEnd, cStart, cEnd, c;

162:   if (!cutLabel) return 0;
163:   DMPlexGetVTKCellHeight(dm, &cellHeight);
164:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
165:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
166:   /* Label vertices that should be duplicated */
167:   DMLabelCreate(PETSC_COMM_SELF, "Cut Vertices", cutVertexLabel);
168:   DMLabelGetStratumIS(cutLabel, 2, &cutcells);
169:   if (cutcells) {
170:     PetscInt n;

172:     ISGetIndices(cutcells, &cutc);
173:     ISGetLocalSize(cutcells, &n);
174:     for (c = 0; c < n; ++c) {
175:       if ((cutc[c] >= cStart) && (cutc[c] < cEnd)) {
176:         PetscInt *closure = NULL;
177:         PetscInt  closureSize, cl, value;

179:         DMPlexGetTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
180:         for (cl = 0; cl < closureSize*2; cl += 2) {
181:           if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) {
182:             DMLabelGetValue(cutLabel, closure[cl], &value);
183:             if (value == 1) {
184:               DMLabelSetValue(*cutVertexLabel, closure[cl], 1);
185:             }
186:           }
187:         }
188:         DMPlexRestoreTransitiveClosure(dm, cutc[c], PETSC_TRUE, &closureSize, &closure);
189:       }
190:     }
191:     ISRestoreIndices(cutcells, &cutc);
192:   }
193:   ISDestroy(&cutcells);
194:   return 0;
195: }

197: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
198: {
199:   DM                      dm;
200:   DM                      dmBC;
201:   PetscSection            section, sectionGlobal;
202:   Vec                     gv;
203:   const char             *name;
204:   PetscViewerVTKFieldType ft;
205:   PetscViewerFormat       format;
206:   PetscInt                seqnum;
207:   PetscReal               seqval;
208:   PetscBool               isseq;

210:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
211:   VecGetDM(v, &dm);
212:   DMGetLocalSection(dm, &section);
213:   DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
214:   DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
215:   if (seqnum >= 0) {
216:     PetscViewerHDF5PushTimestepping(viewer);
217:     PetscViewerHDF5SetTimestep(viewer, seqnum);
218:   }
219:   PetscViewerGetFormat(viewer, &format);
220:   DMGetOutputDM(dm, &dmBC);
221:   DMGetGlobalSection(dmBC, &sectionGlobal);
222:   DMGetGlobalVector(dmBC, &gv);
223:   PetscObjectGetName((PetscObject) v, &name);
224:   PetscObjectSetName((PetscObject) gv, name);
225:   DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
226:   DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
227:   PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
228:   if (isseq) VecView_Seq(gv, viewer);
229:   else       VecView_MPI(gv, viewer);
230:   if (format == PETSC_VIEWER_HDF5_VIZ) {
231:     /* Output visualization representation */
232:     PetscInt numFields, f;
233:     DMLabel  cutLabel, cutVertexLabel = NULL;

235:     PetscSectionGetNumFields(section, &numFields);
236:     DMGetLabel(dm, "periodic_cut", &cutLabel);
237:     for (f = 0; f < numFields; ++f) {
238:       Vec         subv;
239:       IS          is;
240:       const char *fname, *fgroup, *componentName;
241:       char        subname[PETSC_MAX_PATH_LEN];
242:       PetscInt    pStart, pEnd, Nc, c;

244:       DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
245:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
246:       PetscSectionGetFieldName(section, f, &fname);
247:       if (!fname) continue;
248:       PetscViewerHDF5PushGroup(viewer, fgroup);
249:       if (cutLabel) {
250:         const PetscScalar *ga;
251:         PetscScalar       *suba;
252:         PetscInt          gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;

254:         DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
255:         PetscSectionGetFieldComponents(section, f, &Nc);
256:         for (p = pStart; p < pEnd; ++p) {
257:           PetscInt gdof, fdof = 0, val;

259:           PetscSectionGetDof(sectionGlobal, p, &gdof);
260:           if (gdof > 0) PetscSectionGetFieldDof(section, p, f, &fdof);
261:           subSize += fdof;
262:           DMLabelGetValue(cutVertexLabel, p, &val);
263:           if (val == 1) extSize += fdof;
264:         }
265:         VecCreate(PetscObjectComm((PetscObject) gv), &subv);
266:         VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
267:         VecSetBlockSize(subv, Nc);
268:         VecSetType(subv, VECSTANDARD);
269:         VecGetOwnershipRange(gv, &gstart, NULL);
270:         VecGetArrayRead(gv, &ga);
271:         VecGetArray(subv, &suba);
272:         for (p = pStart; p < pEnd; ++p) {
273:           PetscInt gdof, goff, val;

275:           PetscSectionGetDof(sectionGlobal, p, &gdof);
276:           if (gdof > 0) {
277:             PetscInt fdof, fc, f2, poff = 0;

279:             PetscSectionGetOffset(sectionGlobal, p, &goff);
280:             /* Can get rid of this loop by storing field information in the global section */
281:             for (f2 = 0; f2 < f; ++f2) {
282:               PetscSectionGetFieldDof(section, p, f2, &fdof);
283:               poff += fdof;
284:             }
285:             PetscSectionGetFieldDof(section, p, f, &fdof);
286:             for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
287:             DMLabelGetValue(cutVertexLabel, p, &val);
288:             if (val == 1) {
289:               for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
290:             }
291:           }
292:         }
293:         VecRestoreArrayRead(gv, &ga);
294:         VecRestoreArray(subv, &suba);
295:         DMLabelDestroy(&cutVertexLabel);
296:       } else {
297:         PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
298:       }
299:       PetscStrncpy(subname, name,sizeof(subname));
300:       PetscStrlcat(subname, "_",sizeof(subname));
301:       PetscStrlcat(subname, fname,sizeof(subname));
302:       PetscObjectSetName((PetscObject) subv, subname);
303:       if (isseq) VecView_Seq(subv, viewer);
304:       else       VecView_MPI(subv, viewer);
305:       if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
306:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "vector");
307:       } else {
308:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, "vector_field_type", PETSC_STRING, "scalar");
309:       }

311:       /* Output the component names in the field if available */
312:       PetscSectionGetFieldComponents(section, f, &Nc);
313:       for (c = 0; c < Nc; ++c){
314:         char componentNameLabel[PETSC_MAX_PATH_LEN];
315:         PetscSectionGetComponentName(section, f, c, &componentName);
316:         PetscSNPrintf(componentNameLabel, sizeof(componentNameLabel), "componentName%D", c);
317:         PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) subv, componentNameLabel, PETSC_STRING, componentName);
318:       }

320:       if (cutLabel) VecDestroy(&subv);
321:       else          PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
322:       PetscViewerHDF5PopGroup(viewer);
323:     }
324:   }
325:   if (seqnum >= 0) {
326:     PetscViewerHDF5PopTimestepping(viewer);
327:   }
328:   DMRestoreGlobalVector(dmBC, &gv);
329:   return 0;
330: }

332: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
333: {
334:   DM             dm;
335:   Vec            locv;
336:   PetscObject    isZero;
337:   const char    *name;
338:   PetscReal      time;

340:   VecGetDM(v, &dm);
341:   DMGetLocalVector(dm, &locv);
342:   PetscObjectGetName((PetscObject) v, &name);
343:   PetscObjectSetName((PetscObject) locv, name);
344:   PetscObjectQuery((PetscObject) v, "__Vec_bc_zero__", &isZero);
345:   PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", isZero);
346:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
347:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
348:   DMGetOutputSequenceNumber(dm, NULL, &time);
349:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
350:   PetscViewerHDF5PushGroup(viewer, "/fields");
351:   PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
352:   VecView_Plex_Local_HDF5_Internal(locv, viewer);
353:   PetscViewerPopFormat(viewer);
354:   PetscViewerHDF5PopGroup(viewer);
355:   PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", NULL);
356:   DMRestoreLocalVector(dm, &locv);
357:   return 0;
358: }

360: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
361: {
362:   PetscBool      isseq;

364:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
365:   PetscViewerHDF5PushGroup(viewer, "/fields");
366:   if (isseq) VecView_Seq(v, viewer);
367:   else       VecView_MPI(v, viewer);
368:   PetscViewerHDF5PopGroup(viewer);
369:   return 0;
370: }

372: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
373: {
374:   DM             dm;
375:   Vec            locv;
376:   const char    *name;
377:   PetscInt       seqnum;

379:   VecGetDM(v, &dm);
380:   DMGetLocalVector(dm, &locv);
381:   PetscObjectGetName((PetscObject) v, &name);
382:   PetscObjectSetName((PetscObject) locv, name);
383:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
384:   PetscViewerHDF5PushGroup(viewer, "/fields");
385:   if (seqnum >= 0) {
386:     PetscViewerHDF5PushTimestepping(viewer);
387:     PetscViewerHDF5SetTimestep(viewer, seqnum);
388:   }
389:   VecLoad_Plex_Local(locv, viewer);
390:   if (seqnum >= 0) {
391:     PetscViewerHDF5PopTimestepping(viewer);
392:   }
393:   PetscViewerHDF5PopGroup(viewer);
394:   DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
395:   DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
396:   DMRestoreLocalVector(dm, &locv);
397:   return 0;
398: }

400: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
401: {
402:   DM             dm;
403:   PetscInt       seqnum;

405:   VecGetDM(v, &dm);
406:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
407:   PetscViewerHDF5PushGroup(viewer, "/fields");
408:   if (seqnum >= 0) {
409:     PetscViewerHDF5PushTimestepping(viewer);
410:     PetscViewerHDF5SetTimestep(viewer, seqnum);
411:   }
412:   VecLoad_Default(v, viewer);
413:   if (seqnum >= 0) {
414:     PetscViewerHDF5PopTimestepping(viewer);
415:   }
416:   PetscViewerHDF5PopGroup(viewer);
417:   return 0;
418: }

420: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
421: {
422:   const char           *topologydm_name;
423:   const char           *pointsName, *coneSizesName, *conesName, *orientationsName;
424:   IS                    pointsIS, coneSizesIS, conesIS, orientationsIS;
425:   PetscInt             *points, *coneSizes, *cones, *orientations;
426:   const PetscInt       *gpoint;
427:   PetscInt              pStart, pEnd, nPoints = 0, conesSize = 0;
428:   PetscInt              p, c, s;
429:   DMPlexStorageVersion  version;
430:   char                  group[PETSC_MAX_PATH_LEN];
431:   MPI_Comm              comm;

433:   PetscObjectGetComm((PetscObject) dm, &comm);
434:   pointsName        = "order";
435:   coneSizesName     = "cones";
436:   conesName         = "cells";
437:   orientationsName  = "orientation";
438:   DMPlexStorageVersionSetUpWriting_Private(dm, viewer, &version);
439:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
440:   if (version.major <= 1) {
441:     PetscStrcpy(group, "/topology");
442:   } else {
443:     PetscSNPrintf(group, sizeof(group), "topologies/%s/topology", topologydm_name);
444:   }
445:   PetscViewerHDF5PushGroup(viewer, group);
446:   DMPlexGetChart(dm, &pStart, &pEnd);
447:   ISGetIndices(globalPointNumbers, &gpoint);
448:   for (p = pStart; p < pEnd; ++p) {
449:     if (gpoint[p] >= 0) {
450:       PetscInt coneSize;

452:       DMPlexGetConeSize(dm, p, &coneSize);
453:       nPoints += 1;
454:       conesSize += coneSize;
455:     }
456:   }
457:   PetscMalloc1(nPoints, &points);
458:   PetscMalloc1(nPoints, &coneSizes);
459:   PetscMalloc1(conesSize, &cones);
460:   PetscMalloc1(conesSize, &orientations);
461:   for (p = pStart, c = 0, s = 0; p < pEnd; ++p) {
462:     if (gpoint[p] >= 0) {
463:       const PetscInt *cone, *ornt;
464:       PetscInt        coneSize, cp;

466:       DMPlexGetConeSize(dm, p, &coneSize);
467:       DMPlexGetCone(dm, p, &cone);
468:       DMPlexGetConeOrientation(dm, p, &ornt);
469:       points[s]    = gpoint[p];
470:       coneSizes[s] = coneSize;
471:       for (cp = 0; cp < coneSize; ++cp, ++c) {
472:         cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]];
473:         orientations[c] = ornt[cp];
474:       }
475:       ++s;
476:     }
477:   }
480:   ISCreateGeneral(comm, nPoints, points, PETSC_OWN_POINTER, &pointsIS);
481:   ISCreateGeneral(comm, nPoints, coneSizes, PETSC_OWN_POINTER, &coneSizesIS);
482:   ISCreateGeneral(comm, conesSize, cones, PETSC_OWN_POINTER, &conesIS);
483:   ISCreateGeneral(comm, conesSize, orientations, PETSC_OWN_POINTER, &orientationsIS);
484:   PetscObjectSetName((PetscObject) pointsIS, pointsName);
485:   PetscObjectSetName((PetscObject) coneSizesIS, coneSizesName);
486:   PetscObjectSetName((PetscObject) conesIS, conesName);
487:   PetscObjectSetName((PetscObject) orientationsIS, orientationsName);
488:   ISView(pointsIS, viewer);
489:   ISView(coneSizesIS, viewer);
490:   ISView(conesIS, viewer);
491:   ISView(orientationsIS, viewer);
492:   ISDestroy(&pointsIS);
493:   ISDestroy(&coneSizesIS);
494:   ISDestroy(&conesIS);
495:   ISDestroy(&orientationsIS);
496:   ISRestoreIndices(globalPointNumbers, &gpoint);
497:   {
498:     PetscInt dim;

500:     DMGetDimension(dm, &dim);
501:     PetscViewerHDF5WriteAttribute(viewer, conesName, "cell_dim", PETSC_INT, (void *) &dim);
502:   }
503:   PetscViewerHDF5PopGroup(viewer);
504:   return 0;
505: }

507: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
508: {
509:   PetscSF         sfPoint;
510:   DMLabel         cutLabel, cutVertexLabel = NULL;
511:   IS              globalVertexNumbers, cutvertices = NULL;
512:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
513:   PetscInt       *vertices;
514:   PetscInt        conesSize = 0;
515:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;

517:   *numCorners = 0;
518:   DMGetDimension(dm, &dim);
519:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
520:   ISGetIndices(globalCellNumbers, &gcell);

522:   for (cell = cStart; cell < cEnd; ++cell) {
523:     PetscInt *closure = NULL;
524:     PetscInt  closureSize, v, Nc = 0;

526:     if (gcell[cell] < 0) continue;
527:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
528:     for (v = 0; v < closureSize*2; v += 2) {
529:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
530:     }
531:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
532:     conesSize += Nc;
533:     if (!numCornersLocal)           numCornersLocal = Nc;
534:     else if (numCornersLocal != Nc) numCornersLocal = 1;
535:   }
536:   MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
538:   /* Handle periodic cuts by identifying vertices which should be duplicated */
539:   DMGetLabel(dm, "periodic_cut", &cutLabel);
540:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
541:   if (cutVertexLabel) DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);
542:   if (cutvertices) {
543:     ISGetIndices(cutvertices, &cutverts);
544:     ISGetLocalSize(cutvertices, &vExtra);
545:   }
546:   DMGetPointSF(dm, &sfPoint);
547:   if (cutLabel) {
548:     const PetscInt    *ilocal;
549:     const PetscSFNode *iremote;
550:     PetscInt           nroots, nleaves;

552:     PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
553:     if (nleaves < 0) {
554:       PetscObjectReference((PetscObject) sfPoint);
555:     } else {
556:       PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
557:       PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, (PetscInt*)ilocal, PETSC_COPY_VALUES, (PetscSFNode*)iremote, PETSC_COPY_VALUES);
558:     }
559:   } else {
560:     PetscObjectReference((PetscObject) sfPoint);
561:   }
562:   /* Number all vertices */
563:   DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
564:   PetscSFDestroy(&sfPoint);
565:   /* Create cones */
566:   ISGetIndices(globalVertexNumbers, &gvertex);
567:   PetscMalloc1(conesSize, &vertices);
568:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
569:     PetscInt *closure = NULL;
570:     PetscInt  closureSize, Nc = 0, p, value = -1;
571:     PetscBool replace;

573:     if (gcell[cell] < 0) continue;
574:     if (cutLabel) DMLabelGetValue(cutLabel, cell, &value);
575:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
576:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
577:     for (p = 0; p < closureSize*2; p += 2) {
578:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
579:         closure[Nc++] = closure[p];
580:       }
581:     }
582:     DMPlexReorderCell(dm, cell, closure);
583:     for (p = 0; p < Nc; ++p) {
584:       PetscInt nv, gv = gvertex[closure[p] - vStart];

586:       if (replace) {
587:         PetscFindInt(closure[p], vExtra, cutverts, &nv);
588:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
589:       }
590:       vertices[v++] = gv < 0 ? -(gv+1) : gv;
591:     }
592:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
593:   }
594:   ISRestoreIndices(globalVertexNumbers, &gvertex);
595:   ISDestroy(&globalVertexNumbers);
596:   ISRestoreIndices(globalCellNumbers, &gcell);
597:   if (cutvertices) ISRestoreIndices(cutvertices, &cutverts);
598:   ISDestroy(&cutvertices);
599:   DMLabelDestroy(&cutVertexLabel);
601:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
602:   PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
603:   PetscObjectSetName((PetscObject) *cellIS, "cells");
604:   return 0;
605: }

607: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
608: {
609:   DM              cdm;
610:   DMLabel         depthLabel, ctLabel;
611:   IS              cellIS;
612:   PetscInt        dim, depth, cellHeight, c;
613:   hid_t           fileId, groupId;

615:   PetscViewerHDF5PushGroup(viewer, "/viz");
616:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
617:   PetscStackCallHDF5(H5Gclose,(groupId));

619:   PetscViewerHDF5PopGroup(viewer);
620:   DMGetDimension(dm, &dim);
621:   DMPlexGetDepth(dm, &depth);
622:   DMGetCoordinateDM(dm, &cdm);
623:   DMPlexGetVTKCellHeight(dm, &cellHeight);
624:   DMPlexGetDepthLabel(dm, &depthLabel);
625:   DMPlexGetCellTypeLabel(dm, &ctLabel);
626:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
627:     const DMPolytopeType ict = (DMPolytopeType) c;
628:     PetscInt             pStart, pEnd, dep, numCorners, n = 0;
629:     PetscBool            output = PETSC_FALSE, doOutput;

631:     if (ict == DM_POLYTOPE_FV_GHOST) continue;
632:     DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
633:     if (pStart >= 0) {
634:       DMLabelGetValue(depthLabel, pStart, &dep);
635:       if (dep == depth - cellHeight) output = PETSC_TRUE;
636:     }
637:     MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
638:     if (!doOutput) continue;
639:     CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners,  &cellIS);
640:     if (!n) {
641:       PetscViewerHDF5PushGroup(viewer, "/viz/topology");
642:     } else {
643:       char group[PETSC_MAX_PATH_LEN];

645:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
646:       PetscViewerHDF5PushGroup(viewer, group);
647:     }
648:     ISView(cellIS, viewer);
649:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
650:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim",     PETSC_INT, (void *) &dim);
651:     ISDestroy(&cellIS);
652:     PetscViewerHDF5PopGroup(viewer);
653:     ++n;
654:   }
655:   return 0;
656: }

658: static PetscErrorCode DMPlexCoordinatesView_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
659: {
660:   DM             cdm;
661:   Vec            coordinates, newcoords;
662:   PetscReal      lengthScale;
663:   PetscInt       m, M, bs;

665:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
666:   DMGetCoordinateDM(dm, &cdm);
667:   DMGetCoordinates(dm, &coordinates);
668:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
669:   PetscObjectSetName((PetscObject) newcoords, "vertices");
670:   VecGetSize(coordinates, &M);
671:   VecGetLocalSize(coordinates, &m);
672:   VecSetSizes(newcoords, m, M);
673:   VecGetBlockSize(coordinates, &bs);
674:   VecSetBlockSize(newcoords, bs);
675:   VecSetType(newcoords,VECSTANDARD);
676:   VecCopy(coordinates, newcoords);
677:   VecScale(newcoords, lengthScale);
678:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
679:   PetscViewerHDF5PushGroup(viewer, "/geometry");
680:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
681:   VecView(newcoords, viewer);
682:   PetscViewerPopFormat(viewer);
683:   PetscViewerHDF5PopGroup(viewer);
684:   VecDestroy(&newcoords);
685:   return 0;
686: }

688: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
689: {
690:   DM              cdm;
691:   Vec             coords, newcoords;
692:   PetscInt        m, M, bs;
693:   PetscReal       lengthScale;
694:   const char     *topologydm_name, *coordinatedm_name, *coordinates_name;

696:   {
697:     PetscViewerFormat     format;
698:     DMPlexStorageVersion  version;

700:     PetscViewerGetFormat(viewer, &format);
701:     DMPlexStorageVersionSetUpWriting_Private(dm, viewer, &version);
702:     if (format == PETSC_VIEWER_HDF5_XDMF || format == PETSC_VIEWER_HDF5_VIZ || version.major <= 1) {
703:       DMPlexCoordinatesView_HDF5_Legacy_Private(dm, viewer);
704:       return 0;
705:     }
706:   }
707:   DMGetCoordinateDM(dm, &cdm);
708:   DMGetCoordinates(dm, &coords);
709:   PetscObjectGetName((PetscObject)cdm, &coordinatedm_name);
710:   PetscObjectGetName((PetscObject)coords, &coordinates_name);
711:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
712:   PetscViewerHDF5PushGroup(viewer, "topologies");
713:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
714:   PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING, coordinatedm_name);
715:   PetscViewerHDF5WriteAttribute(viewer, NULL, "coordinatesName", PETSC_STRING, coordinates_name);
716:   PetscViewerHDF5PopGroup(viewer);
717:   PetscViewerHDF5PopGroup(viewer);
718:   DMPlexSectionView(dm, viewer, cdm);
719:   VecCreate(PetscObjectComm((PetscObject)coords), &newcoords);
720:   PetscObjectSetName((PetscObject)newcoords, coordinates_name);
721:   VecGetSize(coords, &M);
722:   VecGetLocalSize(coords, &m);
723:   VecSetSizes(newcoords, m, M);
724:   VecGetBlockSize(coords, &bs);
725:   VecSetBlockSize(newcoords, bs);
726:   VecSetType(newcoords,VECSTANDARD);
727:   VecCopy(coords, newcoords);
728:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
729:   VecScale(newcoords, lengthScale);
730:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
731:   DMPlexGlobalVectorView(dm, viewer, cdm, newcoords);
732:   PetscViewerPopFormat(viewer);
733:   VecDestroy(&newcoords);
734:   return 0;
735: }

737: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
738: {
739:   DM               cdm;
740:   Vec              coordinatesLocal, newcoords;
741:   PetscSection     cSection, cGlobalSection;
742:   PetscScalar     *coords, *ncoords;
743:   DMLabel          cutLabel, cutVertexLabel = NULL;
744:   const PetscReal *L;
745:   const DMBoundaryType *bd;
746:   PetscReal        lengthScale;
747:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
748:   PetscBool        localized, embedded;
749:   hid_t            fileId, groupId;

751:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
752:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
753:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
754:   VecGetBlockSize(coordinatesLocal, &bs);
755:   DMGetCoordinatesLocalized(dm, &localized);
756:   if (localized == PETSC_FALSE) return 0;
757:   DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
758:   DMGetCoordinateDM(dm, &cdm);
759:   DMGetLocalSection(cdm, &cSection);
760:   DMGetGlobalSection(cdm, &cGlobalSection);
761:   DMGetLabel(dm, "periodic_cut", &cutLabel);
762:   N    = 0;

764:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
765:   VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
766:   PetscSectionGetDof(cSection, vStart, &dof);
767:   PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
768:   embedded  = (PetscBool) (L && dof == 2 && !cutLabel);
769:   if (cutVertexLabel) {
770:     DMLabelGetStratumSize(cutVertexLabel, 1, &v);
771:     N   += dof*v;
772:   }
773:   for (v = vStart; v < vEnd; ++v) {
774:     PetscSectionGetDof(cGlobalSection, v, &dof);
775:     if (dof < 0) continue;
776:     if (embedded) N += dof+1;
777:     else          N += dof;
778:   }
779:   if (embedded) VecSetBlockSize(newcoords, bs+1);
780:   else          VecSetBlockSize(newcoords, bs);
781:   VecSetSizes(newcoords, N, PETSC_DETERMINE);
782:   VecSetType(newcoords, VECSTANDARD);
783:   VecGetArray(coordinatesLocal, &coords);
784:   VecGetArray(newcoords,        &ncoords);
785:   coordSize = 0;
786:   for (v = vStart; v < vEnd; ++v) {
787:     PetscSectionGetDof(cGlobalSection, v, &dof);
788:     PetscSectionGetOffset(cSection, v, &off);
789:     if (dof < 0) continue;
790:     if (embedded) {
791:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
792:         PetscReal theta, phi, r, R;
793:         /* XY-periodic */
794:         /* Suppose its an y-z circle, then
795:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
796:            and the circle in that plane is
797:              \hat r cos(phi) + \hat x sin(phi) */
798:         theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
799:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
800:         r     = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
801:         R     = L[1]/(2.0*PETSC_PI);
802:         ncoords[coordSize++] =  PetscSinReal(phi) * r;
803:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
804:         ncoords[coordSize++] =  PetscSinReal(theta) * (R + r * PetscCosReal(phi));
805:       } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
806:         /* X-periodic */
807:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
808:         ncoords[coordSize++] = coords[off+1];
809:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
810:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
811:         /* Y-periodic */
812:         ncoords[coordSize++] = coords[off+0];
813:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
814:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
815:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
816:         PetscReal phi, r, R;
817:         /* Mobius strip */
818:         /* Suppose its an x-z circle, then
819:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
820:            and in that plane we rotate by pi as we go around the circle
821:              \hat r cos(phi/2) + \hat y sin(phi/2) */
822:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
823:         R     = L[0];
824:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
825:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
826:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
827:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
828:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
829:     } else {
830:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
831:     }
832:   }
833:   if (cutVertexLabel) {
834:     IS              vertices;
835:     const PetscInt *verts;
836:     PetscInt        n;

838:     DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
839:     if (vertices) {
840:       ISGetIndices(vertices, &verts);
841:       ISGetLocalSize(vertices, &n);
842:       for (v = 0; v < n; ++v) {
843:         PetscSectionGetDof(cSection, verts[v], &dof);
844:         PetscSectionGetOffset(cSection, verts[v], &off);
845:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
846:       }
847:       ISRestoreIndices(vertices, &verts);
848:       ISDestroy(&vertices);
849:     }
850:   }
852:   DMLabelDestroy(&cutVertexLabel);
853:   VecRestoreArray(coordinatesLocal, &coords);
854:   VecRestoreArray(newcoords,        &ncoords);
855:   PetscObjectSetName((PetscObject) newcoords, "vertices");
856:   VecScale(newcoords, lengthScale);
857:   PetscViewerHDF5PushGroup(viewer, "/viz");
858:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
859:   PetscStackCallHDF5(H5Gclose,(groupId));
860:   PetscViewerHDF5PopGroup(viewer);
861:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
862:   VecView(newcoords, viewer);
863:   PetscViewerHDF5PopGroup(viewer);
864:   VecDestroy(&newcoords);
865:   return 0;
866: }

868: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
869: {
870:   const char           *topologydm_name;
871:   const PetscInt       *gpoint;
872:   PetscInt              numLabels, l;
873:   DMPlexStorageVersion  version;
874:   char                  group[PETSC_MAX_PATH_LEN];

876:   DMPlexStorageVersionSetUpWriting_Private(dm, viewer, &version);
877:   ISGetIndices(globalPointNumbers, &gpoint);
878:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
879:   if (version.major <= 1) {
880:     PetscStrcpy(group, "/labels");
881:   } else {
882:     PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
883:   }
884:   PetscViewerHDF5PushGroup(viewer, group);
885:   DMGetNumLabels(dm, &numLabels);
886:   for (l = 0; l < numLabels; ++l) {
887:     DMLabel         label;
888:     const char     *name;
889:     IS              valueIS, pvalueIS, globalValueIS;
890:     const PetscInt *values;
891:     PetscInt        numValues, v;
892:     PetscBool       isDepth, output;

894:     DMGetLabelByNum(dm, l, &label);
895:     PetscObjectGetName((PetscObject)label, &name);
896:     DMGetLabelOutput(dm, name, &output);
897:     PetscStrncmp(name, "depth", 10, &isDepth);
898:     if (isDepth || !output) continue;
899:     PetscViewerHDF5PushGroup(viewer, name);
900:     DMLabelGetValueIS(label, &valueIS);
901:     /* Must copy to a new IS on the global comm */
902:     ISGetLocalSize(valueIS, &numValues);
903:     ISGetIndices(valueIS, &values);
904:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
905:     ISRestoreIndices(valueIS, &values);
906:     ISAllGather(pvalueIS, &globalValueIS);
907:     ISDestroy(&pvalueIS);
908:     ISSortRemoveDups(globalValueIS);
909:     ISGetLocalSize(globalValueIS, &numValues);
910:     ISGetIndices(globalValueIS, &values);
911:     for (v = 0; v < numValues; ++v) {
912:       IS              stratumIS, globalStratumIS;
913:       const PetscInt *spoints = NULL;
914:       PetscInt       *gspoints, n = 0, gn, p;
915:       const char     *iname = "indices";
916:       char            group[PETSC_MAX_PATH_LEN];

918:       PetscSNPrintf(group, sizeof(group), "%D", values[v]);
919:       PetscViewerHDF5PushGroup(viewer, group);
920:       DMLabelGetStratumIS(label, values[v], &stratumIS);

922:       if (stratumIS) ISGetLocalSize(stratumIS, &n);
923:       if (stratumIS) ISGetIndices(stratumIS, &spoints);
924:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
925:       PetscMalloc1(gn,&gspoints);
926:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
927:       if (stratumIS) ISRestoreIndices(stratumIS, &spoints);
928:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
929:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

931:       ISView(globalStratumIS, viewer);
932:       ISDestroy(&globalStratumIS);
933:       ISDestroy(&stratumIS);
934:       PetscViewerHDF5PopGroup(viewer);
935:     }
936:     ISRestoreIndices(globalValueIS, &values);
937:     ISDestroy(&globalValueIS);
938:     ISDestroy(&valueIS);
939:     PetscViewerHDF5PopGroup(viewer);
940:   }
941:   ISRestoreIndices(globalPointNumbers, &gpoint);
942:   PetscViewerHDF5PopGroup(viewer);
943:   return 0;
944: }

946: /* We only write cells and vertices. Does this screw up parallel reading? */
947: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
948: {
949:   IS                globalPointNumbers;
950:   PetscViewerFormat format;
951:   PetscBool         viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;

953:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
954:   DMPlexCoordinatesView_HDF5_Internal(dm, viewer);

956:   PetscViewerGetFormat(viewer, &format);
957:   switch (format) {
958:     case PETSC_VIEWER_HDF5_VIZ:
959:       viz_geom    = PETSC_TRUE;
960:       xdmf_topo   = PETSC_TRUE;
961:       break;
962:     case PETSC_VIEWER_HDF5_XDMF:
963:       xdmf_topo   = PETSC_TRUE;
964:       break;
965:     case PETSC_VIEWER_HDF5_PETSC:
966:       petsc_topo  = PETSC_TRUE;
967:       break;
968:     case PETSC_VIEWER_DEFAULT:
969:     case PETSC_VIEWER_NATIVE:
970:       viz_geom    = PETSC_TRUE;
971:       xdmf_topo   = PETSC_TRUE;
972:       petsc_topo  = PETSC_TRUE;
973:       break;
974:     default:
975:       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
976:   }

978:   if (viz_geom)   DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);
979:   if (xdmf_topo)  DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);
980:   if (petsc_topo) {
981:     DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer);
982:     DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer);
983:   }

985:   ISDestroy(&globalPointNumbers);
986:   return 0;
987: }

989: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
990: {
991:   MPI_Comm       comm;
992:   const char    *topologydm_name;
993:   const char    *sectiondm_name;
994:   PetscSection   gsection;

996:   PetscObjectGetComm((PetscObject)sectiondm, &comm);
997:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
998:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
999:   PetscViewerHDF5PushGroup(viewer, "topologies");
1000:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1001:   PetscViewerHDF5PushGroup(viewer, "dms");
1002:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1003:   DMGetGlobalSection(sectiondm, &gsection);
1004:   /* Save raw section */
1005:   PetscSectionView(gsection, viewer);
1006:   /* Save plex wrapper */
1007:   {
1008:     PetscInt        pStart, pEnd, p, n;
1009:     IS              globalPointNumbers;
1010:     const PetscInt *gpoints;
1011:     IS              orderIS;
1012:     PetscInt       *order;

1014:     PetscSectionGetChart(gsection, &pStart, &pEnd);
1015:     DMPlexCreatePointNumbering(dm, &globalPointNumbers);
1016:     ISGetIndices(globalPointNumbers, &gpoints);
1017:     for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) n++;
1018:     /* "order" is an array of global point numbers.
1019:        When loading, it is used with topology/order array
1020:        to match section points with plex topology points. */
1021:     PetscMalloc1(n, &order);
1022:     for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) order[n++] = gpoints[p];
1023:     ISRestoreIndices(globalPointNumbers, &gpoints);
1024:     ISDestroy(&globalPointNumbers);
1025:     ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS);
1026:     PetscObjectSetName((PetscObject)orderIS, "order");
1027:     ISView(orderIS, viewer);
1028:     ISDestroy(&orderIS);
1029:   }
1030:   PetscViewerHDF5PopGroup(viewer);
1031:   PetscViewerHDF5PopGroup(viewer);
1032:   PetscViewerHDF5PopGroup(viewer);
1033:   PetscViewerHDF5PopGroup(viewer);
1034:   return 0;
1035: }

1037: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1038: {
1039:   const char     *topologydm_name;
1040:   const char     *sectiondm_name;
1041:   const char     *vec_name;
1042:   PetscInt        bs;

1044:   /* Check consistency */
1045:   {
1046:     PetscSF   pointsf, pointsf1;

1048:     DMGetPointSF(dm, &pointsf);
1049:     DMGetPointSF(sectiondm, &pointsf1);
1051:   }
1052:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1053:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1054:   PetscObjectGetName((PetscObject)vec, &vec_name);
1055:   PetscViewerHDF5PushGroup(viewer, "topologies");
1056:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1057:   PetscViewerHDF5PushGroup(viewer, "dms");
1058:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1059:   PetscViewerHDF5PushGroup(viewer, "vecs");
1060:   PetscViewerHDF5PushGroup(viewer, vec_name);
1061:   VecGetBlockSize(vec, &bs);
1062:   PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
1063:   VecSetBlockSize(vec, 1);
1064:   /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but,    */
1065:   /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view    */
1066:   /* is set to VecView_Plex, which would save vec in a predefined location.  */
1067:   /* To save vec in where we want, we create a new Vec (temp) with           */
1068:   /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
1069:   {
1070:     Vec                temp;
1071:     const PetscScalar *array;
1072:     PetscLayout        map;

1074:     VecCreate(PetscObjectComm((PetscObject)vec), &temp);
1075:     PetscObjectSetName((PetscObject)temp, vec_name);
1076:     VecGetLayout(vec, &map);
1077:     VecSetLayout(temp, map);
1078:     VecSetUp(temp);
1079:     VecGetArrayRead(vec, &array);
1080:     VecPlaceArray(temp, array);
1081:     VecView(temp, viewer);
1082:     VecResetArray(temp);
1083:     VecRestoreArrayRead(vec, &array);
1084:     VecDestroy(&temp);
1085:   }
1086:   VecSetBlockSize(vec, bs);
1087:   PetscViewerHDF5PopGroup(viewer);
1088:   PetscViewerHDF5PopGroup(viewer);
1089:   PetscViewerHDF5PopGroup(viewer);
1090:   PetscViewerHDF5PopGroup(viewer);
1091:   PetscViewerHDF5PopGroup(viewer);
1092:   PetscViewerHDF5PopGroup(viewer);
1093:   return 0;
1094: }

1096: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
1097: {
1098:   MPI_Comm        comm;
1099:   const char     *topologydm_name;
1100:   const char     *sectiondm_name;
1101:   const char     *vec_name;
1102:   PetscSection    section;
1103:   PetscBool       includesConstraints;
1104:   Vec             gvec;
1105:   PetscInt        m, bs;

1107:   PetscObjectGetComm((PetscObject)dm, &comm);
1108:   /* Check consistency */
1109:   {
1110:     PetscSF   pointsf, pointsf1;

1112:     DMGetPointSF(dm, &pointsf);
1113:     DMGetPointSF(sectiondm, &pointsf1);
1115:   }
1116:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1117:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1118:   PetscObjectGetName((PetscObject)vec, &vec_name);
1119:   PetscViewerHDF5PushGroup(viewer, "topologies");
1120:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1121:   PetscViewerHDF5PushGroup(viewer, "dms");
1122:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1123:   PetscViewerHDF5PushGroup(viewer, "vecs");
1124:   PetscViewerHDF5PushGroup(viewer, vec_name);
1125:   VecGetBlockSize(vec, &bs);
1126:   PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
1127:   VecCreate(comm, &gvec);
1128:   PetscObjectSetName((PetscObject)gvec, vec_name);
1129:   DMGetGlobalSection(sectiondm, &section);
1130:   PetscSectionGetIncludesConstraints(section, &includesConstraints);
1131:   if (includesConstraints) PetscSectionGetStorageSize(section, &m);
1132:   else PetscSectionGetConstrainedStorageSize(section, &m);
1133:   VecSetSizes(gvec, m, PETSC_DECIDE);
1134:   VecSetUp(gvec);
1135:   DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec);
1136:   DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec);
1137:   VecView(gvec, viewer);
1138:   VecDestroy(&gvec);
1139:   PetscViewerHDF5PopGroup(viewer);
1140:   PetscViewerHDF5PopGroup(viewer);
1141:   PetscViewerHDF5PopGroup(viewer);
1142:   PetscViewerHDF5PopGroup(viewer);
1143:   PetscViewerHDF5PopGroup(viewer);
1144:   PetscViewerHDF5PopGroup(viewer);
1145:   return 0;
1146: }

1148: struct _n_LoadLabelsCtx {
1149:   MPI_Comm    comm;
1150:   PetscMPIInt rank;
1151:   DM          dm;
1152:   PetscViewer viewer;
1153:   DMLabel     label;
1154:   PetscSF     sfXC;
1155:   PetscLayout layoutX;
1156: };
1157: typedef struct _n_LoadLabelsCtx *LoadLabelsCtx;

1159: static PetscErrorCode LoadLabelsCtxCreate(DM dm, PetscViewer viewer, PetscSF sfXC, LoadLabelsCtx *ctx)
1160: {
1161:   PetscNew(ctx);
1162:   PetscObjectReference((PetscObject) ((*ctx)->dm = dm));
1163:   PetscObjectReference((PetscObject) ((*ctx)->viewer = viewer));
1164:   PetscObjectGetComm((PetscObject)dm, &(*ctx)->comm);
1165:   MPI_Comm_rank((*ctx)->comm, &(*ctx)->rank);
1166:   (*ctx)->sfXC = sfXC;
1167:   if (sfXC) {
1168:     PetscInt nX;

1170:     PetscObjectReference((PetscObject) sfXC);
1171:     PetscSFGetGraph(sfXC, &nX, NULL, NULL, NULL);
1172:     PetscLayoutCreateFromSizes((*ctx)->comm, nX, PETSC_DECIDE, 1, &(*ctx)->layoutX);
1173:   }
1174:   return 0;
1175: }

1177: static PetscErrorCode LoadLabelsCtxDestroy(LoadLabelsCtx *ctx)
1178: {
1179:   if (!*ctx) return 0;
1180:   DMDestroy(&(*ctx)->dm);
1181:   PetscViewerDestroy(&(*ctx)->viewer);
1182:   PetscSFDestroy(&(*ctx)->sfXC);
1183:   PetscLayoutDestroy(&(*ctx)->layoutX);
1184:   PetscFree(*ctx);
1185:   return 0;
1186: }

1188: /*
1189:     A: on-disk points
1190:     X: global points [0, NX)
1191:     C: distributed plex points
1192: */
1193: static herr_t ReadLabelStratumHDF5_Distribute_Private(IS stratumIS, LoadLabelsCtx ctx, IS *newStratumIS)
1194: {
1195:   MPI_Comm        comm    = ctx->comm;
1196:   PetscSF         sfXC    = ctx->sfXC;
1197:   PetscLayout     layoutX = ctx->layoutX;
1198:   PetscSF         sfXA;
1199:   const PetscInt *A_points;
1200:   PetscInt        nX, nC;
1201:   PetscInt        n;

1203:   PetscSFGetGraph(sfXC, &nX, &nC, NULL, NULL);
1204:   ISGetLocalSize(stratumIS, &n);
1205:   ISGetIndices(stratumIS, &A_points);
1206:   PetscSFCreate(comm, &sfXA);
1207:   PetscSFSetGraphLayout(sfXA, layoutX, n, NULL, PETSC_USE_POINTER, A_points);
1208:   ISCreate(comm, newStratumIS);
1209:   ISSetType(*newStratumIS,ISGENERAL);
1210:   {
1211:     PetscInt    i;
1212:     PetscBool  *A_mask, *X_mask, *C_mask;

1214:     PetscCalloc3(n, &A_mask, nX, &X_mask, nC, &C_mask);
1215:     for (i=0; i<n; i++) A_mask[i] = PETSC_TRUE;
1216:     PetscSFReduceBegin(sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1217:     PetscSFReduceEnd(  sfXA, MPIU_BOOL, A_mask, X_mask, MPI_REPLACE);
1218:     PetscSFBcastBegin( sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1219:     PetscSFBcastEnd(   sfXC, MPIU_BOOL, X_mask, C_mask, MPI_LOR);
1220:     ISGeneralSetIndicesFromMask(*newStratumIS, 0, nC, C_mask);
1221:     PetscFree3(A_mask, X_mask, C_mask);
1222:   }
1223:   PetscSFDestroy(&sfXA);
1224:   ISRestoreIndices(stratumIS, &A_points);
1225:   return 0;
1226: }

1228: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *vname, const H5L_info_t *info, void *op_data)
1229: {
1230:   LoadLabelsCtx   ctx    = (LoadLabelsCtx) op_data;
1231:   PetscViewer     viewer = ctx->viewer;
1232:   DMLabel         label  = ctx->label;
1233:   MPI_Comm        comm   = ctx->comm;
1234:   IS              stratumIS;
1235:   const PetscInt *ind;
1236:   PetscInt        value, N, i;

1238:   PetscOptionsStringToInt(vname, &value);
1239:   ISCreate(comm, &stratumIS);
1240:   PetscObjectSetName((PetscObject) stratumIS, "indices");
1241:   PetscViewerHDF5PushGroup(viewer, vname); /* labels/<lname>/<vname> */

1243:   if (!ctx->sfXC) {
1244:     /* Force serial load */
1245:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
1246:     PetscLayoutSetLocalSize(stratumIS->map, !ctx->rank ? N : 0);
1247:     PetscLayoutSetSize(stratumIS->map, N);
1248:   }
1249:   ISLoad(stratumIS, viewer);

1251:   if (ctx->sfXC) {
1252:     IS newStratumIS;

1254:     ReadLabelStratumHDF5_Distribute_Private(stratumIS, ctx, &newStratumIS);
1255:     ISDestroy(&stratumIS);
1256:     stratumIS = newStratumIS;
1257:   }

1259:   PetscViewerHDF5PopGroup(viewer);
1260:   ISGetLocalSize(stratumIS, &N);
1261:   ISGetIndices(stratumIS, &ind);
1262:   for (i = 0; i < N; ++i) DMLabelSetValue(label, ind[i], value);
1263:   ISRestoreIndices(stratumIS, &ind);
1264:   ISDestroy(&stratumIS);
1265:   return 0;
1266: }

1268: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *lname, const H5L_info_t *info, void *op_data)
1269: {
1270:   LoadLabelsCtx  ctx = (LoadLabelsCtx) op_data;
1271:   DM             dm  = ctx->dm;
1272:   hsize_t        idx = 0;
1274:   PetscBool      flg;
1275:   herr_t         err;

1277:   DMHasLabel(dm, lname, &flg);
1278:   if (flg) DMRemoveLabel(dm, lname, NULL);
1279:   DMCreateLabel(dm, lname); if (ierr) return (herr_t) ierr;
1280:   DMGetLabel(dm, lname, &ctx->label); if (ierr) return (herr_t) ierr;
1281:   PetscViewerHDF5PushGroup(ctx->viewer, lname); /* labels/<lname> */
1282:   /* Iterate over the label's strata */
1283:   PetscStackCallHDF5Return(err, H5Literate_by_name, (g_id, lname, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1284:   PetscViewerHDF5PopGroup(ctx->viewer);
1285:   return err;
1286: }

1288: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1289: {
1290:   const char           *topologydm_name;
1291:   LoadLabelsCtx         ctx;
1292:   hsize_t               idx = 0;
1293:   char                  group[PETSC_MAX_PATH_LEN];
1294:   DMPlexStorageVersion  version;
1295:   PetscBool             distributed, hasGroup;

1297:   DMPlexIsDistributed(dm, &distributed);
1298:   if (distributed) {
1300:   }
1301:   LoadLabelsCtxCreate(dm, viewer, sfXC, &ctx);
1302:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1303:   DMPlexStorageVersionGet_Private(dm, viewer, &version);
1304:   if (version.major <= 1) {
1305:     PetscStrcpy(group, "labels");
1306:   } else {
1307:     PetscSNPrintf(group, sizeof(group), "topologies/%s/labels", topologydm_name);
1308:   }
1309:   PetscViewerHDF5PushGroup(viewer, group);
1310:   PetscViewerHDF5HasGroup(viewer, NULL, &hasGroup);
1311:   if (hasGroup) {
1312:     hid_t fileId, groupId;

1314:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
1315:     /* Iterate over labels */
1316:     PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, ctx));
1317:     PetscStackCallHDF5(H5Gclose,(groupId));
1318:   }
1319:   PetscViewerHDF5PopGroup(viewer);
1320:   LoadLabelsCtxDestroy(&ctx);
1321:   return 0;
1322: }

1324: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sf)
1325: {
1326:   MPI_Comm              comm;
1327:   const char           *topologydm_name;
1328:   const char           *pointsName, *coneSizesName, *conesName, *orientationsName;
1329:   IS                    pointsIS, coneSizesIS, conesIS, orientationsIS;
1330:   const PetscInt       *points, *coneSizes, *cones, *orientations;
1331:   PetscInt             *cone, *ornt;
1332:   PetscInt              dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1333:   PetscMPIInt           size, rank;
1334:   char                  group[PETSC_MAX_PATH_LEN];
1335:   DMPlexStorageVersion  version;

1337:   pointsName        = "order";
1338:   coneSizesName     = "cones";
1339:   conesName         = "cells";
1340:   orientationsName  = "orientation";
1341:   PetscObjectGetComm((PetscObject)dm, &comm);
1342:   MPI_Comm_size(comm, &size);
1343:   MPI_Comm_rank(comm, &rank);
1344:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1345:   DMPlexStorageVersionGet_Private(dm, viewer, &version);
1346:   if (version.major <= 1) {
1347:     PetscStrcpy(group, "/topology");
1348:   } else {
1349:     PetscSNPrintf(group, sizeof(group), "topologies/%s/topology", topologydm_name);
1350:   }
1351:   PetscViewerHDF5PushGroup(viewer, group);
1352:   ISCreate(comm, &pointsIS);
1353:   PetscObjectSetName((PetscObject) pointsIS, pointsName);
1354:   ISCreate(comm, &coneSizesIS);
1355:   PetscObjectSetName((PetscObject) coneSizesIS, coneSizesName);
1356:   ISCreate(comm, &conesIS);
1357:   PetscObjectSetName((PetscObject) conesIS, conesName);
1358:   ISCreate(comm, &orientationsIS);
1359:   PetscObjectSetName((PetscObject) orientationsIS, orientationsName);
1360:   PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) conesIS, "cell_dim", PETSC_INT, NULL, &dim);
1361:   DMSetDimension(dm, dim);
1362:   {
1363:     /* Force serial load */
1364:     PetscViewerHDF5ReadSizes(viewer, pointsName, NULL, &Np);
1365:     PetscLayoutSetLocalSize(pointsIS->map, rank == 0 ? Np : 0);
1366:     PetscLayoutSetSize(pointsIS->map, Np);
1367:     pEnd = rank == 0 ? Np : 0;
1368:     PetscViewerHDF5ReadSizes(viewer, coneSizesName, NULL, &Np);
1369:     PetscLayoutSetLocalSize(coneSizesIS->map, rank == 0 ? Np : 0);
1370:     PetscLayoutSetSize(coneSizesIS->map, Np);
1371:     PetscViewerHDF5ReadSizes(viewer, conesName, NULL, &N);
1372:     PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? N : 0);
1373:     PetscLayoutSetSize(conesIS->map, N);
1374:     PetscViewerHDF5ReadSizes(viewer, orientationsName, NULL, &N);
1375:     PetscLayoutSetLocalSize(orientationsIS->map, rank == 0 ? N : 0);
1376:     PetscLayoutSetSize(orientationsIS->map, N);
1377:   }
1378:   ISLoad(pointsIS, viewer);
1379:   ISLoad(coneSizesIS, viewer);
1380:   ISLoad(conesIS, viewer);
1381:   ISLoad(orientationsIS, viewer);
1382:   PetscViewerHDF5PopGroup(viewer);
1383:   /* Create Plex */
1384:   DMPlexSetChart(dm, 0, pEnd);
1385:   ISGetIndices(pointsIS, &points);
1386:   ISGetIndices(coneSizesIS, &coneSizes);
1387:   for (p = 0; p < pEnd; ++p) {
1388:     DMPlexSetConeSize(dm, points[p], coneSizes[p]);
1389:     maxConeSize = PetscMax(maxConeSize, coneSizes[p]);
1390:   }
1391:   DMSetUp(dm);
1392:   ISGetIndices(conesIS, &cones);
1393:   ISGetIndices(orientationsIS, &orientations);
1394:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
1395:   for (p = 0, q = 0; p < pEnd; ++p) {
1396:     for (c = 0; c < coneSizes[p]; ++c, ++q) {
1397:       cone[c] = cones[q];
1398:       ornt[c] = orientations[q];
1399:     }
1400:     DMPlexSetCone(dm, points[p], cone);
1401:     DMPlexSetConeOrientation(dm, points[p], ornt);
1402:   }
1403:   PetscFree2(cone,ornt);
1404:   /* Create global section migration SF */
1405:   if (sf) {
1406:     PetscLayout  layout;
1407:     PetscInt    *globalIndices;

1409:     PetscMalloc1(pEnd, &globalIndices);
1410:     /* plex point == globalPointNumber in this case */
1411:     for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1412:     PetscLayoutCreate(comm, &layout);
1413:     PetscLayoutSetSize(layout, Np);
1414:     PetscLayoutSetBlockSize(layout, 1);
1415:     PetscLayoutSetUp(layout);
1416:     PetscSFCreate(comm, sf);
1417:     PetscSFSetFromOptions(*sf);
1418:     PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices);
1419:     PetscLayoutDestroy(&layout);
1420:     PetscFree(globalIndices);
1421:   }
1422:   /* Clean-up */
1423:   ISRestoreIndices(pointsIS, &points);
1424:   ISRestoreIndices(coneSizesIS, &coneSizes);
1425:   ISRestoreIndices(conesIS, &cones);
1426:   ISRestoreIndices(orientationsIS, &orientations);
1427:   ISDestroy(&pointsIS);
1428:   ISDestroy(&coneSizesIS);
1429:   ISDestroy(&conesIS);
1430:   ISDestroy(&orientationsIS);
1431:   /* Fill in the rest of the topology structure */
1432:   DMPlexSymmetrize(dm);
1433:   DMPlexStratify(dm);
1434:   return 0;
1435: }

1437: /* If the file is old, it not only has different path to the coordinates, but   */
1438: /* does not contain coordinateDMs, so must fall back to the old implementation. */
1439: static PetscErrorCode DMPlexCoordinatesLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1440: {
1441:   PetscSection    coordSection;
1442:   Vec             coordinates;
1443:   PetscReal       lengthScale;
1444:   PetscInt        spatialDim, N, numVertices, vStart, vEnd, v;
1445:   PetscMPIInt     rank;

1447:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1448:   /* Read geometry */
1449:   PetscViewerHDF5PushGroup(viewer, "/geometry");
1450:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
1451:   PetscObjectSetName((PetscObject) coordinates, "vertices");
1452:   {
1453:     /* Force serial load */
1454:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
1455:     VecSetSizes(coordinates, !rank ? N : 0, N);
1456:     VecSetBlockSize(coordinates, spatialDim);
1457:   }
1458:   VecLoad(coordinates, viewer);
1459:   PetscViewerHDF5PopGroup(viewer);
1460:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1461:   VecScale(coordinates, 1.0/lengthScale);
1462:   VecGetLocalSize(coordinates, &numVertices);
1463:   VecGetBlockSize(coordinates, &spatialDim);
1464:   numVertices /= spatialDim;
1465:   /* Create coordinates */
1466:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1468:   DMGetCoordinateSection(dm, &coordSection);
1469:   PetscSectionSetNumFields(coordSection, 1);
1470:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
1471:   PetscSectionSetChart(coordSection, vStart, vEnd);
1472:   for (v = vStart; v < vEnd; ++v) {
1473:     PetscSectionSetDof(coordSection, v, spatialDim);
1474:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1475:   }
1476:   PetscSectionSetUp(coordSection);
1477:   DMSetCoordinates(dm, coordinates);
1478:   VecDestroy(&coordinates);
1479:   return 0;
1480: }

1482: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF sfXC)
1483: {
1484:   DM                    cdm;
1485:   Vec                   coords;
1486:   PetscInt              blockSize;
1487:   PetscReal             lengthScale;
1488:   PetscSF               lsf;
1489:   const char           *topologydm_name;
1490:   char                 *coordinatedm_name, *coordinates_name;

1492:   {
1493:     DMPlexStorageVersion  version;

1495:     DMPlexStorageVersionGet_Private(dm, viewer, &version);
1496:     if (version.major <= 1) {
1497:       DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1498:       return 0;
1499:     }
1500:   }
1502:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1503:   PetscViewerHDF5PushGroup(viewer, "topologies");
1504:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1505:   PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinateDMName", PETSC_STRING , NULL, &coordinatedm_name);
1506:   PetscViewerHDF5ReadAttribute(viewer, NULL, "coordinatesName", PETSC_STRING , NULL, &coordinates_name);
1507:   PetscViewerHDF5PopGroup(viewer);
1508:   PetscViewerHDF5PopGroup(viewer);
1509:   DMGetCoordinateDM(dm, &cdm);
1510:   PetscObjectSetName((PetscObject)cdm, coordinatedm_name);
1511:   PetscFree(coordinatedm_name);
1512:   /* lsf: on-disk data -> in-memory local vector associated with cdm's local section */
1513:   DMPlexSectionLoad(dm, viewer, cdm, sfXC, NULL, &lsf);
1514:   DMCreateLocalVector(cdm, &coords);
1515:   PetscObjectSetName((PetscObject)coords, coordinates_name);
1516:   PetscFree(coordinates_name);
1517:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
1518:   DMPlexLocalVectorLoad(dm, viewer, cdm, lsf, coords);
1519:   PetscViewerPopFormat(viewer);
1520:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1521:   VecScale(coords, 1.0/lengthScale);
1522:   DMSetCoordinatesLocal(dm, coords);
1523:   VecGetBlockSize(coords, &blockSize);
1524:   DMSetCoordinateDim(dm, blockSize);
1525:   VecDestroy(&coords);
1526:   PetscSFDestroy(&lsf);
1527:   return 0;
1528: }

1530: static PetscErrorCode DMPlexLoad_HDF5_Legacy_Private(DM dm, PetscViewer viewer)
1531: {
1532:   DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL);
1533:   DMPlexLabelsLoad_HDF5_Internal(dm, viewer, NULL);
1534:   DMPlexCoordinatesLoad_HDF5_Legacy_Private(dm, viewer);
1535:   return 0;
1536: }

1538: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1539: {
1540:   PetscSF               sfXC;

1542:   {
1543:     DMPlexStorageVersion  version;

1545:     DMPlexStorageVersionGet_Private(dm, viewer, &version);
1546:     if (version.major <= 1) {
1547:       DMPlexLoad_HDF5_Legacy_Private(dm, viewer);
1548:       return 0;
1549:     }
1550:   }
1551:   DMPlexTopologyLoad_HDF5_Internal(dm, viewer, &sfXC);
1552:   DMPlexLabelsLoad_HDF5_Internal(dm, viewer, sfXC);
1553:   DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer, sfXC);
1554:   PetscSFDestroy(&sfXC);
1555:   return 0;
1556: }

1558: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1559: {
1560:   MPI_Comm        comm;
1561:   PetscInt        pStart, pEnd, p, m;
1562:   PetscInt       *goffs, *ilocal;
1563:   PetscBool       rootIncludeConstraints, leafIncludeConstraints;

1565:   PetscObjectGetComm((PetscObject)leafSection, &comm);
1566:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1567:   PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints);
1568:   PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints);
1569:   if (rootIncludeConstraints && leafIncludeConstraints) PetscSectionGetStorageSize(leafSection, &m);
1570:   else PetscSectionGetConstrainedStorageSize(leafSection, &m);
1571:   PetscMalloc1(m, &ilocal);
1572:   PetscMalloc1(m, &goffs);
1573:   /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
1574:   /* for the top-level section (not for each field), so one must have   */
1575:   /* rootSection->pointMajor == PETSC_TRUE.                             */
1577:   /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
1579:   for (p = pStart, m = 0; p < pEnd; ++p) {
1580:     PetscInt        dof, cdof, i, j, off, goff;
1581:     const PetscInt *cinds;

1583:     PetscSectionGetDof(leafSection, p, &dof);
1584:     if (dof < 0) continue;
1585:     goff = globalOffsets[p-pStart];
1586:     PetscSectionGetOffset(leafSection, p, &off);
1587:     PetscSectionGetConstraintDof(leafSection, p, &cdof);
1588:     PetscSectionGetConstraintIndices(leafSection, p, &cinds);
1589:     for (i = 0, j = 0; i < dof; ++i) {
1590:       PetscBool constrained = (PetscBool) (j < cdof && i == cinds[j]);

1592:       if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {ilocal[m] = off++; goffs[m++] = goff++;}
1593:       else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
1594:       else if (!leafIncludeConstraints &&  rootIncludeConstraints) ++goff;
1595:       if (constrained) ++j;
1596:     }
1597:   }
1598:   PetscSFCreate(comm, sectionSF);
1599:   PetscSFSetFromOptions(*sectionSF);
1600:   PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs);
1601:   PetscFree(goffs);
1602:   return 0;
1603: }

1605: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
1606: {
1607:   MPI_Comm       comm;
1608:   PetscMPIInt    size, rank;
1609:   const char    *topologydm_name;
1610:   const char    *sectiondm_name;
1611:   PetscSection   sectionA, sectionB;
1612:   PetscInt       nX, n, i;
1613:   PetscSF        sfAB;

1615:   PetscObjectGetComm((PetscObject)dm, &comm);
1616:   MPI_Comm_size(comm, &size);
1617:   MPI_Comm_rank(comm, &rank);
1618:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1619:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1620:   PetscViewerHDF5PushGroup(viewer, "topologies");
1621:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1622:   PetscViewerHDF5PushGroup(viewer, "dms");
1623:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1624:   /* A: on-disk points                        */
1625:   /* X: list of global point numbers, [0, NX) */
1626:   /* B: plex points                           */
1627:   /* Load raw section (sectionA)              */
1628:   PetscSectionCreate(comm, &sectionA);
1629:   PetscSectionLoad(sectionA, viewer);
1630:   PetscSectionGetChart(sectionA, NULL, &n);
1631:   /* Create sfAB: A -> B */
1632: #if defined(PETSC_USE_DEBUG)
1633:   {
1634:     PetscInt  N, N1;

1636:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1);
1637:     MPI_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm);
1639:   }
1640: #endif
1641:   {
1642:     IS              orderIS;
1643:     const PetscInt *gpoints;
1644:     PetscSF         sfXA, sfAX;
1645:     PetscLayout     layout;
1646:     PetscSFNode    *owners, *buffer;
1647:     PetscInt        nleaves;
1648:     PetscInt       *ilocal;
1649:     PetscSFNode    *iremote;

1651:     /* Create sfAX: A -> X */
1652:     ISCreate(comm, &orderIS);
1653:     PetscObjectSetName((PetscObject)orderIS, "order");
1654:     PetscLayoutSetLocalSize(orderIS->map, n);
1655:     ISLoad(orderIS, viewer);
1656:     PetscLayoutCreate(comm, &layout);
1657:     PetscSFGetGraph(sfXB, &nX, NULL, NULL, NULL);
1658:     PetscLayoutSetLocalSize(layout, nX);
1659:     PetscLayoutSetBlockSize(layout, 1);
1660:     PetscLayoutSetUp(layout);
1661:     PetscSFCreate(comm, &sfXA);
1662:     ISGetIndices(orderIS, &gpoints);
1663:     PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints);
1664:     ISRestoreIndices(orderIS, &gpoints);
1665:     ISDestroy(&orderIS);
1666:     PetscLayoutDestroy(&layout);
1667:     PetscMalloc1(n, &owners);
1668:     PetscMalloc1(nX, &buffer);
1669:     for (i = 0; i < n; ++i) {owners[i].rank = rank; owners[i].index = i;}
1670:     for (i = 0; i < nX; ++i) {buffer[i].rank = -1; buffer[i].index = -1;}
1671:     PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1672:     PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1673:     PetscSFDestroy(&sfXA);
1674:     PetscFree(owners);
1675:     for (i = 0, nleaves = 0; i < nX; ++i) if (buffer[i].rank >= 0) nleaves++;
1676:     PetscMalloc1(nleaves, &ilocal);
1677:     PetscMalloc1(nleaves, &iremote);
1678:     for (i = 0, nleaves = 0; i < nX; ++i) {
1679:       if (buffer[i].rank >= 0) {
1680:         ilocal[nleaves] = i;
1681:         iremote[nleaves].rank = buffer[i].rank;
1682:         iremote[nleaves].index = buffer[i].index;
1683:         nleaves++;
1684:       }
1685:     }
1686:     PetscSFCreate(comm, &sfAX);
1687:     PetscSFSetFromOptions(sfAX);
1688:     PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1689:     /* Fix PetscSFCompose() and replace the code-block below with:  */
1690:     /* PetscSFCompose(sfAX, sfXB, &sfAB);      */
1691:     /* which currently causes segmentation fault due to sparse map. */
1692:     {
1693:       PetscInt     npoints;
1694:       PetscInt     mleaves;
1695:       PetscInt    *jlocal;
1696:       PetscSFNode *jremote;

1698:       PetscSFGetGraph(sfXB, NULL, &npoints, NULL, NULL);
1699:       PetscMalloc1(npoints, &owners);
1700:       for (i = 0; i < npoints; ++i) {owners[i].rank = -1; owners[i].index = -1;}
1701:       PetscSFBcastBegin(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1702:       PetscSFBcastEnd(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1703:       for (i = 0, mleaves = 0; i < npoints; ++i) if (owners[i].rank >= 0) mleaves++;
1704:       PetscMalloc1(mleaves, &jlocal);
1705:       PetscMalloc1(mleaves, &jremote);
1706:       for (i = 0, mleaves = 0; i < npoints; ++i) {
1707:         if (owners[i].rank >= 0) {
1708:           jlocal[mleaves] = i;
1709:           jremote[mleaves].rank = owners[i].rank;
1710:           jremote[mleaves].index = owners[i].index;
1711:           mleaves++;
1712:         }
1713:       }
1714:       PetscSFCreate(comm, &sfAB);
1715:       PetscSFSetFromOptions(sfAB);
1716:       PetscSFSetGraph(sfAB, n, mleaves, jlocal, PETSC_OWN_POINTER, jremote, PETSC_OWN_POINTER);
1717:       PetscFree(owners);
1718:     }
1719:     PetscFree(buffer);
1720:     PetscSFDestroy(&sfAX);
1721:   }
1722:   PetscViewerHDF5PopGroup(viewer);
1723:   PetscViewerHDF5PopGroup(viewer);
1724:   PetscViewerHDF5PopGroup(viewer);
1725:   PetscViewerHDF5PopGroup(viewer);
1726:   /* Create plex section (sectionB) */
1727:   DMGetLocalSection(sectiondm, &sectionB);
1728:   if (lsf || gsf) {
1729:     PetscLayout  layout;
1730:     PetscInt     M, m;
1731:     PetscInt    *offsetsA;
1732:     PetscBool    includesConstraintsA;

1734:     PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB);
1735:     PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA);
1736:     if (includesConstraintsA) PetscSectionGetStorageSize(sectionA, &m);
1737:     else PetscSectionGetConstrainedStorageSize(sectionA, &m);
1738:     MPI_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm);
1739:     PetscLayoutCreate(comm, &layout);
1740:     PetscLayoutSetSize(layout, M);
1741:     PetscLayoutSetUp(layout);
1742:     if (lsf) {
1743:       PetscSF lsfABdata;

1745:       DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata);
1746:       *lsf = lsfABdata;
1747:     }
1748:     if (gsf) {
1749:       PetscSection  gsectionB, gsectionB1;
1750:       PetscBool     includesConstraintsB;
1751:       PetscSF       gsfABdata, pointsf;

1753:       DMGetGlobalSection(sectiondm, &gsectionB1);
1754:       PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB);
1755:       DMGetPointSF(sectiondm, &pointsf);
1756:       PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB);
1757:       DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata);
1758:       PetscSectionDestroy(&gsectionB);
1759:       *gsf = gsfABdata;
1760:     }
1761:     PetscLayoutDestroy(&layout);
1762:     PetscFree(offsetsA);
1763:   } else {
1764:     PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB);
1765:   }
1766:   PetscSFDestroy(&sfAB);
1767:   PetscSectionDestroy(&sectionA);
1768:   return 0;
1769: }

1771: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
1772: {
1773:   MPI_Comm           comm;
1774:   const char        *topologydm_name;
1775:   const char        *sectiondm_name;
1776:   const char        *vec_name;
1777:   Vec                vecA;
1778:   PetscInt           mA, m, bs;
1779:   const PetscInt    *ilocal;
1780:   const PetscScalar *src;
1781:   PetscScalar       *dest;

1783:   PetscObjectGetComm((PetscObject)dm, &comm);
1784:   DMPlexGetHDF5Name_Private(dm, &topologydm_name);
1785:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1786:   PetscObjectGetName((PetscObject)vec, &vec_name);
1787:   PetscViewerHDF5PushGroup(viewer, "topologies");
1788:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1789:   PetscViewerHDF5PushGroup(viewer, "dms");
1790:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1791:   PetscViewerHDF5PushGroup(viewer, "vecs");
1792:   PetscViewerHDF5PushGroup(viewer, vec_name);
1793:   VecCreate(comm, &vecA);
1794:   PetscObjectSetName((PetscObject)vecA, vec_name);
1795:   PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL);
1796:   /* Check consistency */
1797:   {
1798:     PetscSF   pointsf, pointsf1;
1799:     PetscInt  m1, i, j;

1801:     DMGetPointSF(dm, &pointsf);
1802:     DMGetPointSF(sectiondm, &pointsf1);
1804: #if defined(PETSC_USE_DEBUG)
1805:     {
1806:       PetscInt  MA, MA1;

1808:       MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm);
1809:       PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1);
1811:     }
1812: #endif
1813:     VecGetLocalSize(vec, &m1);
1815:     for (i = 0; i < m; ++i) {
1816:       j = ilocal ? ilocal[i] : i;
1818:     }
1819:   }
1820:   VecSetSizes(vecA, mA, PETSC_DECIDE);
1821:   VecLoad(vecA, viewer);
1822:   VecGetArrayRead(vecA, &src);
1823:   VecGetArray(vec, &dest);
1824:   PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1825:   PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1826:   VecRestoreArray(vec, &dest);
1827:   VecRestoreArrayRead(vecA, &src);
1828:   VecDestroy(&vecA);
1829:   PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *) &bs);
1830:   VecSetBlockSize(vec, bs);
1831:   PetscViewerHDF5PopGroup(viewer);
1832:   PetscViewerHDF5PopGroup(viewer);
1833:   PetscViewerHDF5PopGroup(viewer);
1834:   PetscViewerHDF5PopGroup(viewer);
1835:   PetscViewerHDF5PopGroup(viewer);
1836:   PetscViewerHDF5PopGroup(viewer);
1837:   return 0;
1838: }
1839: #endif