Actual source code: plexhdf5.c

petsc-3.11.2 2019-05-18
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petsc/private/isimpl.h>
  3:  #include <petsc/private/vecimpl.h>
  4: #include <petscviewerhdf5.h>

  6: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

260: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
261: {
262:   PetscBool      isseq;

266:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
267:   PetscViewerHDF5PushGroup(viewer, "/fields");
268:   if (isseq) {VecView_Seq(v, viewer);}
269:   else       {VecView_MPI(v, viewer);}
270:   PetscViewerHDF5PopGroup(viewer);
271:   return(0);
272: }

274: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
275: {
276:   DM             dm;
277:   Vec            locv;
278:   const char    *name;
279:   PetscInt       seqnum;

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

298: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
299: {
300:   DM             dm;
301:   PetscInt       seqnum;

305:   VecGetDM(v, &dm);
306:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
307:   PetscViewerHDF5SetTimestep(viewer, seqnum);
308:   PetscViewerHDF5PushGroup(viewer, "/fields");
309:   VecLoad_Default(v, viewer);
310:   PetscViewerHDF5PopGroup(viewer);
311:   return(0);
312: }

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

323:   ISGetIndices(globalPointNumbers, &gpoint);
324:   DMGetDimension(dm, &dim);
325:   DMPlexGetChart(dm, &pStart, &pEnd);
326:   for (p = pStart; p < pEnd; ++p) {
327:     if (gpoint[p] >= 0) {
328:       PetscInt coneSize;

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

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

377: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
378: {
379:   DM              cdm;
380:   PetscSF         sfPoint;
381:   DMLabel         cutLabel, cutVertexLabel = NULL;
382:   PetscSection    cSection;
383:   IS              cellIS, globalVertexNumbers, cutvertices = NULL;
384:   const PetscInt *gvertex, *cutverts = NULL;
385:   PetscInt       *vertices;
386:   PetscInt        dim, depth, vStart, vEnd, vExtra = 0, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
387:   hid_t           fileId, groupId;
388:   PetscErrorCode  ierr;

391:   DMGetDimension(dm, &dim);
392:   DMPlexGetDepth(dm, &depth);
393:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
394:   DMGetCoordinateDM(dm, &cdm);
395:   DMGetDefaultGlobalSection(cdm, &cSection);
396:   DMPlexGetVTKCellHeight(dm, &cellHeight);
397:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
398:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
399:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
400:   for (cell = cStart; cell < cEnd; ++cell) {
401:     PetscInt *closure = NULL;
402:     PetscInt  closureSize, v, Nc = 0;

404:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
405:     for (v = 0; v < closureSize*2; v += 2) {
406:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
407:     }
408:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
409:     conesSize += Nc;
410:     if (!numCornersLocal)           numCornersLocal = Nc;
411:     else if (numCornersLocal != Nc) numCornersLocal = 1;
412:   }
413:   MPIU_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
414:   if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");

416:   PetscViewerHDF5PushGroup(viewer, "/viz");
417:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
418:   PetscStackCallHDF5(H5Gclose,(groupId));
419:   PetscViewerHDF5PopGroup(viewer);

421:   DMGetLabel(dm, "periodic_cut", &cutLabel);
422:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
423:   if (cutVertexLabel) {DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);}
424:   if (cutvertices) {
425:     ISGetIndices(cutvertices, &cutverts);
426:     ISGetLocalSize(cutvertices, &vExtra);
427:   }

429:   DMGetPointSF(dm, &sfPoint);
430:   if (cutLabel) {
431:     const PetscInt    *ilocal;
432:     const PetscSFNode *iremote;
433:     PetscInt           nroots, nleaves;

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

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

462:       if (replace) {
463:         PetscFindInt(closure[p], vExtra, cutverts, &nv);
464:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
465:       }
466:       vertices[v++] = gv < 0 ? -(gv+1) : gv;
467:     }
468:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
469:   }
470:   ISRestoreIndices(globalVertexNumbers, &gvertex);
471:   ISDestroy(&globalVertexNumbers);
472:   if (cutvertices) {ISRestoreIndices(cutvertices, &cutverts);}
473:   ISDestroy(&cutvertices);
474:   DMLabelDestroy(&cutVertexLabel);
475:   if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
476:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);
477:   PetscLayoutSetBlockSize(cellIS->map, numCorners);
478:   PetscObjectSetName((PetscObject) cellIS, "cells");
479:   PetscViewerHDF5PushGroup(viewer, "/viz/topology");
480:   ISView(cellIS, viewer);
481:   PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
482:   PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim",     PETSC_INT, (void *) &dim);
483:   ISDestroy(&cellIS);
484:   PetscViewerHDF5PopGroup(viewer);
485:   return(0);
486: }

488: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
489: {
490:   DM             cdm;
491:   Vec            coordinates, newcoords;
492:   PetscReal      lengthScale;
493:   PetscInt       m, M, bs;

497:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
498:   DMGetCoordinateDM(dm, &cdm);
499:   DMGetCoordinates(dm, &coordinates);
500:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
501:   PetscObjectSetName((PetscObject) newcoords, "vertices");
502:   VecGetSize(coordinates, &M);
503:   VecGetLocalSize(coordinates, &m);
504:   VecSetSizes(newcoords, m, M);
505:   VecGetBlockSize(coordinates, &bs);
506:   VecSetBlockSize(newcoords, bs);
507:   VecSetType(newcoords,VECSTANDARD);
508:   VecCopy(coordinates, newcoords);
509:   VecScale(newcoords, lengthScale);
510:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
511:   PetscViewerHDF5PushGroup(viewer, "/geometry");
512:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
513:   VecView(newcoords, viewer);
514:   PetscViewerPopFormat(viewer);
515:   PetscViewerHDF5PopGroup(viewer);
516:   VecDestroy(&newcoords);
517:   return(0);
518: }

520: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
521: {
522:   DM               cdm;
523:   Vec              coordinatesLocal, newcoords;
524:   PetscSection     cSection, cGlobalSection;
525:   PetscScalar     *coords, *ncoords;
526:   DMLabel          cutLabel, cutVertexLabel = NULL;
527:   const PetscReal *L;
528:   const DMBoundaryType *bd;
529:   PetscReal        lengthScale;
530:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
531:   PetscBool        localized, embedded;
532:   hid_t            fileId, groupId;
533:   PetscErrorCode   ierr;

536:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
537:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
538:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
539:   VecGetBlockSize(coordinatesLocal, &bs);
540:   DMGetCoordinatesLocalized(dm, &localized);
541:   if (localized == PETSC_FALSE) return(0);
542:   DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
543:   DMGetCoordinateDM(dm, &cdm);
544:   DMGetDefaultSection(cdm, &cSection);
545:   DMGetDefaultGlobalSection(cdm, &cGlobalSection);
546:   DMGetLabel(dm, "periodic_cut", &cutLabel);
547:   N    = 0;

549:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
550:   VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
551:   PetscSectionGetDof(cSection, vStart, &dof);
552:   PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
553:   embedded  = (PetscBool) (L && dof == 2 && !cutLabel);
554:   if (cutVertexLabel) {
555:     DMLabelGetStratumSize(cutVertexLabel, 1, &v);
556:     N   += dof*v;
557:   }
558:   for (v = vStart; v < vEnd; ++v) {
559:     PetscSectionGetDof(cGlobalSection, v, &dof);
560:     if (dof < 0) continue;
561:     if (embedded) N += dof+1;
562:     else          N += dof;
563:   }
564:   if (embedded) {VecSetBlockSize(newcoords, bs+1);}
565:   else          {VecSetBlockSize(newcoords, bs);}
566:   VecSetSizes(newcoords, N, PETSC_DETERMINE);
567:   VecSetType(newcoords, VECSTANDARD);
568:   VecGetArray(coordinatesLocal, &coords);
569:   VecGetArray(newcoords,        &ncoords);
570:   coordSize = 0;
571:   for (v = vStart; v < vEnd; ++v) {
572:     PetscSectionGetDof(cGlobalSection, v, &dof);
573:     PetscSectionGetOffset(cSection, v, &off);
574:     if (dof < 0) continue;
575:     if (embedded) {
576:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
577:         PetscReal theta, phi, r, R;
578:         /* XY-periodic */
579:         /* Suppose its an y-z circle, then
580:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
581:            and the circle in that plane is
582:              \hat r cos(phi) + \hat x sin(phi) */
583:         theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
584:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
585:         r     = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
586:         R     = L[1]/(2.0*PETSC_PI);
587:         ncoords[coordSize++] =  PetscSinReal(phi) * r;
588:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
589:         ncoords[coordSize++] =  PetscSinReal(theta) * (R + r * PetscCosReal(phi));
590:       } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
591:         /* X-periodic */
592:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
593:         ncoords[coordSize++] = coords[off+1];
594:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
595:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
596:         /* Y-periodic */
597:         ncoords[coordSize++] = coords[off+0];
598:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
599:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
600:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
601:         PetscReal phi, r, R;
602:         /* Mobius strip */
603:         /* Suppose its an x-z circle, then
604:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
605:            and in that plane we rotate by pi as we go around the circle
606:              \hat r cos(phi/2) + \hat y sin(phi/2) */
607:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
608:         R     = L[0];
609:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
610:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
611:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
612:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
613:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
614:     } else {
615:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
616:     }
617:   }
618:   if (cutVertexLabel) {
619:     IS              vertices;
620:     const PetscInt *verts;
621:     PetscInt        n;

623:     DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
624:     if (vertices) {
625:       ISGetIndices(vertices, &verts);
626:       ISGetLocalSize(vertices, &n);
627:       for (v = 0; v < n; ++v) {
628:         PetscSectionGetDof(cSection, verts[v], &dof);
629:         PetscSectionGetOffset(cSection, verts[v], &off);
630:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
631:       }
632:       ISRestoreIndices(vertices, &verts);
633:       ISDestroy(&vertices);
634:     }
635:   }
636:   if (coordSize != N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %D != %D", coordSize, N);
637:   DMLabelDestroy(&cutVertexLabel);
638:   VecRestoreArray(coordinatesLocal, &coords);
639:   VecRestoreArray(newcoords,        &ncoords);
640:   PetscObjectSetName((PetscObject) newcoords, "vertices");
641:   VecScale(newcoords, lengthScale);
642:   PetscViewerHDF5PushGroup(viewer, "/viz");
643:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
644:   PetscStackCallHDF5(H5Gclose,(groupId));
645:   PetscViewerHDF5PopGroup(viewer);
646:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
647:   VecView(newcoords, viewer);
648:   PetscViewerHDF5PopGroup(viewer);
649:   VecDestroy(&newcoords);
650:   return(0);
651: }


654: static PetscErrorCode DMPlexWriteLabels_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
655: {
656:   const PetscInt   *gpoint;
657:   PetscInt          numLabels, l;
658:   hid_t             fileId, groupId;
659:   PetscErrorCode    ierr;

662:   ISGetIndices(globalPointNumbers, &gpoint);
663:   PetscViewerHDF5PushGroup(viewer, "/labels");
664:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
665:   if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
666:   PetscViewerHDF5PopGroup(viewer);
667:   DMGetNumLabels(dm, &numLabels);
668:   for (l = 0; l < numLabels; ++l) {
669:     DMLabel         label;
670:     const char     *name;
671:     IS              valueIS, pvalueIS, globalValueIS;
672:     const PetscInt *values;
673:     PetscInt        numValues, v;
674:     PetscBool       isDepth, output;
675:     char            group[PETSC_MAX_PATH_LEN];

677:     DMGetLabelName(dm, l, &name);
678:     DMGetLabelOutput(dm, name, &output);
679:     PetscStrncmp(name, "depth", 10, &isDepth);
680:     if (isDepth || !output) continue;
681:     DMGetLabel(dm, name, &label);
682:     PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
683:     PetscViewerHDF5PushGroup(viewer, group);
684:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
685:     if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
686:     PetscViewerHDF5PopGroup(viewer);
687:     DMLabelGetValueIS(label, &valueIS);
688:     /* Must copy to a new IS on the global comm */
689:     ISGetLocalSize(valueIS, &numValues);
690:     ISGetIndices(valueIS, &values);
691:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
692:     ISRestoreIndices(valueIS, &values);
693:     ISAllGather(pvalueIS, &globalValueIS);
694:     ISDestroy(&pvalueIS);
695:     ISSortRemoveDups(globalValueIS);
696:     ISGetLocalSize(globalValueIS, &numValues);
697:     ISGetIndices(globalValueIS, &values);
698:     for (v = 0; v < numValues; ++v) {
699:       IS              stratumIS, globalStratumIS;
700:       const PetscInt *spoints = NULL;
701:       PetscInt       *gspoints, n = 0, gn, p;
702:       const char     *iname = "indices";

704:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%d", name, values[v]);
705:       DMLabelGetStratumIS(label, values[v], &stratumIS);

707:       if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
708:       if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
709:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
710:       PetscMalloc1(gn,&gspoints);
711:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
712:       if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
713:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
714:       if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
715:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

717:       PetscViewerHDF5PushGroup(viewer, group);
718:       ISView(globalStratumIS, viewer);
719:       PetscViewerHDF5PopGroup(viewer);
720:       ISDestroy(&globalStratumIS);
721:       ISDestroy(&stratumIS);
722:     }
723:     ISRestoreIndices(globalValueIS, &values);
724:     ISDestroy(&globalValueIS);
725:     ISDestroy(&valueIS);
726:   }
727:   ISRestoreIndices(globalPointNumbers, &gpoint);
728:   return(0);
729: }

731: /* We only write cells and vertices. Does this screw up parallel reading? */
732: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
733: {
734:   IS                globalPointNumbers;
735:   PetscViewerFormat format;
736:   PetscBool         viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
737:   PetscErrorCode    ierr;

740:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
741:   DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
742:   DMPlexWriteLabels_HDF5_Static(dm, globalPointNumbers, viewer);

744:   PetscViewerGetFormat(viewer, &format);
745:   switch (format) {
746:     case PETSC_VIEWER_HDF5_VIZ:
747:       viz_geom    = PETSC_TRUE;
748:       xdmf_topo   = PETSC_TRUE;
749:       break;
750:     case PETSC_VIEWER_HDF5_XDMF:
751:       xdmf_topo   = PETSC_TRUE;
752:       break;
753:     case PETSC_VIEWER_HDF5_PETSC:
754:       petsc_topo  = PETSC_TRUE;
755:       break;
756:     case PETSC_VIEWER_DEFAULT:
757:     case PETSC_VIEWER_NATIVE:
758:       viz_geom    = PETSC_TRUE;
759:       xdmf_topo   = PETSC_TRUE;
760:       petsc_topo  = PETSC_TRUE;
761:       break;
762:     default:
763:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
764:   }

766:   if (viz_geom)   {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
767:   if (xdmf_topo)  {DMPlexWriteTopology_Vertices_HDF5_Static(dm, viewer);}
768:   if (petsc_topo) {DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);}

770:   ISDestroy(&globalPointNumbers);
771:   return(0);
772: }

774: typedef struct {
775:   PetscMPIInt rank;
776:   DM          dm;
777:   PetscViewer viewer;
778:   DMLabel     label;
779: } LabelCtx;

781: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
782: {
783:   PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
784:   DMLabel         label  = ((LabelCtx *) op_data)->label;
785:   IS              stratumIS;
786:   const PetscInt *ind;
787:   PetscInt        value, N, i;
788:   const char     *lname;
789:   char            group[PETSC_MAX_PATH_LEN];
790:   PetscErrorCode  ierr;

792:   PetscOptionsStringToInt(name, &value);
793:   ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
794:   PetscObjectSetName((PetscObject) stratumIS, "indices");
795:   PetscObjectGetName((PetscObject) label, &lname);
796:   PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
797:   PetscViewerHDF5PushGroup(viewer, group);
798:   {
799:     /* Force serial load */
800:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
801:     PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
802:     PetscLayoutSetSize(stratumIS->map, N);
803:   }
804:   ISLoad(stratumIS, viewer);
805:   PetscViewerHDF5PopGroup(viewer);
806:   ISGetLocalSize(stratumIS, &N);
807:   ISGetIndices(stratumIS, &ind);
808:   for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
809:   ISRestoreIndices(stratumIS, &ind);
810:   ISDestroy(&stratumIS);
811:   return 0;
812: }

814: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
815: {
816:   DM             dm  = ((LabelCtx *) op_data)->dm;
817:   hsize_t        idx = 0;
819:   herr_t         err;

821:   DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
822:   DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
823:   PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
824:   return err;
825: }

827: PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM dm, PetscViewer viewer)
828: {
829:   LabelCtx        ctx;
830:   hid_t           fileId, groupId;
831:   hsize_t         idx = 0;
832:   PetscErrorCode  ierr;

835:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
836:   ctx.dm     = dm;
837:   ctx.viewer = viewer;
838:   PetscViewerHDF5PushGroup(viewer, "/labels");
839:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
840:   PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
841:   PetscStackCallHDF5(H5Gclose,(groupId));
842:   PetscViewerHDF5PopGroup(viewer);
843:   return(0);
844: }

846: /* The first version will read everything onto proc 0, letting the user distribute
847:    The next will create a naive partition, and then rebalance after reading
848: */
849: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
850: {
851:   PetscSection    coordSection;
852:   Vec             coordinates;
853:   IS              orderIS, conesIS, cellsIS, orntsIS;
854:   const PetscInt *order, *cones, *cells, *ornts;
855:   PetscReal       lengthScale;
856:   PetscInt       *cone, *ornt;
857:   PetscInt        dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
858:   PetscMPIInt     rank;
859:   PetscErrorCode  ierr;

862:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
863:   /* Read toplogy */
864:   PetscViewerHDF5PushGroup(viewer, "/topology");
865:   ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
866:   PetscObjectSetName((PetscObject) orderIS, "order");
867:   ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
868:   PetscObjectSetName((PetscObject) conesIS, "cones");
869:   ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
870:   PetscObjectSetName((PetscObject) cellsIS, "cells");
871:   ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
872:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
873:   PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
874:   DMSetDimension(dm, dim);
875:   {
876:     /* Force serial load */
877:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
878:     PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
879:     PetscLayoutSetSize(orderIS->map, pEnd);
880:     PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
881:     PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
882:     PetscLayoutSetSize(conesIS->map, pEnd);
883:     pEnd = !rank ? pEnd : 0;
884:     PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
885:     PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
886:     PetscLayoutSetSize(cellsIS->map, N);
887:     PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
888:     PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
889:     PetscLayoutSetSize(orntsIS->map, N);
890:   }
891:   ISLoad(orderIS, viewer);
892:   ISLoad(conesIS, viewer);
893:   ISLoad(cellsIS, viewer);
894:   ISLoad(orntsIS, viewer);
895:   PetscViewerHDF5PopGroup(viewer);
896:   /* Read geometry */
897:   PetscViewerHDF5PushGroup(viewer, "/geometry");
898:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
899:   PetscObjectSetName((PetscObject) coordinates, "vertices");
900:   {
901:     /* Force serial load */
902:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
903:     VecSetSizes(coordinates, !rank ? N : 0, N);
904:     VecSetBlockSize(coordinates, spatialDim);
905:   }
906:   VecLoad(coordinates, viewer);
907:   PetscViewerHDF5PopGroup(viewer);
908:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
909:   VecScale(coordinates, 1.0/lengthScale);
910:   VecGetLocalSize(coordinates, &numVertices);
911:   VecGetBlockSize(coordinates, &spatialDim);
912:   numVertices /= spatialDim;
913:   /* Create Plex */
914:   DMPlexSetChart(dm, 0, pEnd);
915:   ISGetIndices(orderIS, &order);
916:   ISGetIndices(conesIS, &cones);
917:   for (p = 0; p < pEnd; ++p) {
918:     DMPlexSetConeSize(dm, order[p], cones[p]);
919:     maxConeSize = PetscMax(maxConeSize, cones[p]);
920:   }
921:   DMSetUp(dm);
922:   ISGetIndices(cellsIS, &cells);
923:   ISGetIndices(orntsIS, &ornts);
924:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
925:   for (p = 0, q = 0; p < pEnd; ++p) {
926:     for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
927:     DMPlexSetCone(dm, order[p], cone);
928:     DMPlexSetConeOrientation(dm, order[p], ornt);
929:   }
930:   PetscFree2(cone,ornt);
931:   ISRestoreIndices(orderIS, &order);
932:   ISRestoreIndices(conesIS, &cones);
933:   ISRestoreIndices(cellsIS, &cells);
934:   ISRestoreIndices(orntsIS, &ornts);
935:   ISDestroy(&orderIS);
936:   ISDestroy(&conesIS);
937:   ISDestroy(&cellsIS);
938:   ISDestroy(&orntsIS);
939:   DMPlexSymmetrize(dm);
940:   DMPlexStratify(dm);
941:   /* Create coordinates */
942:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
943:   if (numVertices != vEnd - vStart) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of coordinates loaded %d does not match number of vertices %d", numVertices, vEnd - vStart);
944:   DMGetCoordinateSection(dm, &coordSection);
945:   PetscSectionSetNumFields(coordSection, 1);
946:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
947:   PetscSectionSetChart(coordSection, vStart, vEnd);
948:   for (v = vStart; v < vEnd; ++v) {
949:     PetscSectionSetDof(coordSection, v, spatialDim);
950:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
951:   }
952:   PetscSectionSetUp(coordSection);
953:   DMSetCoordinates(dm, coordinates);
954:   VecDestroy(&coordinates);
955:   /* Read Labels */
956:   DMPlexLoadLabels_HDF5_Internal(dm, viewer);
957:   return(0);
958: }
959: #endif