Actual source code: plexhdf5.c

petsc-3.14.1 2020-11-03
Report Typos and Errors
  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: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

265: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
266: {
267:   PetscBool      isseq;

271:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
272:   PetscViewerHDF5PushGroup(viewer, "/fields");
273:   if (isseq) {VecView_Seq(v, viewer);}
274:   else       {VecView_MPI(v, viewer);}
275:   PetscViewerHDF5PopGroup(viewer);
276:   return(0);
277: }

279: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
280: {
281:   DM             dm;
282:   Vec            locv;
283:   const char    *name;
284:   PetscInt       seqnum;

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

303: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
304: {
305:   DM             dm;
306:   PetscInt       seqnum;

310:   VecGetDM(v, &dm);
311:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
312:   PetscViewerHDF5SetTimestep(viewer, seqnum);
313:   PetscViewerHDF5PushGroup(viewer, "/fields");
314:   VecLoad_Default(v, viewer);
315:   PetscViewerHDF5PopGroup(viewer);
316:   return(0);
317: }

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

328:   ISGetIndices(globalPointNumbers, &gpoint);
329:   DMGetDimension(dm, &dim);
330:   DMPlexGetChart(dm, &pStart, &pEnd);
331:   for (p = pStart; p < pEnd; ++p) {
332:     if (gpoint[p] >= 0) {
333:       PetscInt coneSize;

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

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

382: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
383: {
384:   PetscSF         sfPoint;
385:   DMLabel         cutLabel, cutVertexLabel = NULL;
386:   IS              globalVertexNumbers, cutvertices = NULL;
387:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
388:   PetscInt       *vertices;
389:   PetscInt        conesSize = 0;
390:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
391:   PetscErrorCode  ierr;

394:   *numCorners = 0;
395:   DMGetDimension(dm, &dim);
396:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
397:   ISGetIndices(globalCellNumbers, &gcell);

399:   for (cell = cStart; cell < cEnd; ++cell) {
400:     PetscInt *closure = NULL;
401:     PetscInt  closureSize, v, Nc = 0;

403:     if (gcell[cell] < 0) continue;
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");
415:   /* Handle periodic cuts by identifying vertices which should be duplicated */
416:   DMGetLabel(dm, "periodic_cut", &cutLabel);
417:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
418:   if (cutVertexLabel) {DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);}
419:   if (cutvertices) {
420:     ISGetIndices(cutvertices, &cutverts);
421:     ISGetLocalSize(cutvertices, &vExtra);
422:   }
423:   DMGetPointSF(dm, &sfPoint);
424:   if (cutLabel) {
425:     const PetscInt    *ilocal;
426:     const PetscSFNode *iremote;
427:     PetscInt           nroots, nleaves;

429:     PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
430:     if (nleaves < 0) {
431:       PetscObjectReference((PetscObject) sfPoint);
432:     } else {
433:       PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
434:       PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, ilocal, PETSC_USE_POINTER, iremote, PETSC_USE_POINTER);
435:     }
436:   } else {
437:     PetscObjectReference((PetscObject) sfPoint);
438:   }
439:   /* Number all vertices */
440:   DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
441:   PetscSFDestroy(&sfPoint);
442:   /* Create cones */
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 (gcell[cell] < 0) continue;
451:     if (cutLabel) {DMLabelGetValue(cutLabel, cell, &value);}
452:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
453:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
454:     for (p = 0; p < closureSize*2; p += 2) {
455:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
456:         closure[Nc++] = closure[p];
457:       }
458:     }
459:     DMPlexReorderCell(dm, cell, closure);
460:     for (p = 0; p < Nc; ++p) {
461:       PetscInt nv, gv = gvertex[closure[p] - vStart];

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

484: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
485: {
486:   DM              cdm;
487:   DMLabel         depthLabel, ctLabel;
488:   IS              cellIS;
489:   PetscInt        dim, depth, cellHeight, c;
490:   hid_t           fileId, groupId;
491:   PetscErrorCode  ierr;

494:   PetscViewerHDF5PushGroup(viewer, "/viz");
495:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
496:   PetscStackCallHDF5(H5Gclose,(groupId));

498:   PetscViewerHDF5PopGroup(viewer);
499:   DMGetDimension(dm, &dim);
500:   DMPlexGetDepth(dm, &depth);
501:   DMGetCoordinateDM(dm, &cdm);
502:   DMPlexGetVTKCellHeight(dm, &cellHeight);
503:   DMPlexGetDepthLabel(dm, &depthLabel);
504:   DMPlexGetCellTypeLabel(dm, &ctLabel);
505:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
506:     const DMPolytopeType ict = (DMPolytopeType) c;
507:     PetscInt             pStart, pEnd, dep, numCorners, n = 0;
508:     PetscBool            output = PETSC_FALSE, doOutput;

510:     DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
511:     if (pStart >= 0) {
512:       DMLabelGetValue(depthLabel, pStart, &dep);
513:       if (dep == depth - cellHeight) output = PETSC_TRUE;
514:     }
515:     MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
516:     if (!doOutput) continue;
517:     CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners,  &cellIS);
518:     if (!n) {
519:       PetscViewerHDF5PushGroup(viewer, "/viz/topology");
520:     } else {
521:       char group[PETSC_MAX_PATH_LEN];

523:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
524:       PetscViewerHDF5PushGroup(viewer, group);
525:     }
526:     ISView(cellIS, viewer);
527:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
528:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim",     PETSC_INT, (void *) &dim);
529:     ISDestroy(&cellIS);
530:     PetscViewerHDF5PopGroup(viewer);
531:     ++n;
532:   }
533:   return(0);
534: }

536: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
537: {
538:   DM             cdm;
539:   Vec            coordinates, newcoords;
540:   PetscReal      lengthScale;
541:   PetscInt       m, M, bs;

545:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
546:   DMGetCoordinateDM(dm, &cdm);
547:   DMGetCoordinates(dm, &coordinates);
548:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
549:   PetscObjectSetName((PetscObject) newcoords, "vertices");
550:   VecGetSize(coordinates, &M);
551:   VecGetLocalSize(coordinates, &m);
552:   VecSetSizes(newcoords, m, M);
553:   VecGetBlockSize(coordinates, &bs);
554:   VecSetBlockSize(newcoords, bs);
555:   VecSetType(newcoords,VECSTANDARD);
556:   VecCopy(coordinates, newcoords);
557:   VecScale(newcoords, lengthScale);
558:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
559:   PetscViewerHDF5PushGroup(viewer, "/geometry");
560:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
561:   VecView(newcoords, viewer);
562:   PetscViewerPopFormat(viewer);
563:   PetscViewerHDF5PopGroup(viewer);
564:   VecDestroy(&newcoords);
565:   return(0);
566: }

568: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
569: {
570:   DM               cdm;
571:   Vec              coordinatesLocal, newcoords;
572:   PetscSection     cSection, cGlobalSection;
573:   PetscScalar     *coords, *ncoords;
574:   DMLabel          cutLabel, cutVertexLabel = NULL;
575:   const PetscReal *L;
576:   const DMBoundaryType *bd;
577:   PetscReal        lengthScale;
578:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
579:   PetscBool        localized, embedded;
580:   hid_t            fileId, groupId;
581:   PetscErrorCode   ierr;

584:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
585:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
586:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
587:   VecGetBlockSize(coordinatesLocal, &bs);
588:   DMGetCoordinatesLocalized(dm, &localized);
589:   if (localized == PETSC_FALSE) return(0);
590:   DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
591:   DMGetCoordinateDM(dm, &cdm);
592:   DMGetLocalSection(cdm, &cSection);
593:   DMGetGlobalSection(cdm, &cGlobalSection);
594:   DMGetLabel(dm, "periodic_cut", &cutLabel);
595:   N    = 0;

597:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
598:   VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
599:   PetscSectionGetDof(cSection, vStart, &dof);
600:   PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
601:   embedded  = (PetscBool) (L && dof == 2 && !cutLabel);
602:   if (cutVertexLabel) {
603:     DMLabelGetStratumSize(cutVertexLabel, 1, &v);
604:     N   += dof*v;
605:   }
606:   for (v = vStart; v < vEnd; ++v) {
607:     PetscSectionGetDof(cGlobalSection, v, &dof);
608:     if (dof < 0) continue;
609:     if (embedded) N += dof+1;
610:     else          N += dof;
611:   }
612:   if (embedded) {VecSetBlockSize(newcoords, bs+1);}
613:   else          {VecSetBlockSize(newcoords, bs);}
614:   VecSetSizes(newcoords, N, PETSC_DETERMINE);
615:   VecSetType(newcoords, VECSTANDARD);
616:   VecGetArray(coordinatesLocal, &coords);
617:   VecGetArray(newcoords,        &ncoords);
618:   coordSize = 0;
619:   for (v = vStart; v < vEnd; ++v) {
620:     PetscSectionGetDof(cGlobalSection, v, &dof);
621:     PetscSectionGetOffset(cSection, v, &off);
622:     if (dof < 0) continue;
623:     if (embedded) {
624:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
625:         PetscReal theta, phi, r, R;
626:         /* XY-periodic */
627:         /* Suppose its an y-z circle, then
628:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
629:            and the circle in that plane is
630:              \hat r cos(phi) + \hat x sin(phi) */
631:         theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
632:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
633:         r     = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
634:         R     = L[1]/(2.0*PETSC_PI);
635:         ncoords[coordSize++] =  PetscSinReal(phi) * r;
636:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
637:         ncoords[coordSize++] =  PetscSinReal(theta) * (R + r * PetscCosReal(phi));
638:       } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
639:         /* X-periodic */
640:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
641:         ncoords[coordSize++] = coords[off+1];
642:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
643:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
644:         /* Y-periodic */
645:         ncoords[coordSize++] = coords[off+0];
646:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
647:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
648:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
649:         PetscReal phi, r, R;
650:         /* Mobius strip */
651:         /* Suppose its an x-z circle, then
652:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
653:            and in that plane we rotate by pi as we go around the circle
654:              \hat r cos(phi/2) + \hat y sin(phi/2) */
655:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
656:         R     = L[0];
657:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
658:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
659:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
660:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
661:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
662:     } else {
663:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
664:     }
665:   }
666:   if (cutVertexLabel) {
667:     IS              vertices;
668:     const PetscInt *verts;
669:     PetscInt        n;

671:     DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
672:     if (vertices) {
673:       ISGetIndices(vertices, &verts);
674:       ISGetLocalSize(vertices, &n);
675:       for (v = 0; v < n; ++v) {
676:         PetscSectionGetDof(cSection, verts[v], &dof);
677:         PetscSectionGetOffset(cSection, verts[v], &off);
678:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
679:       }
680:       ISRestoreIndices(vertices, &verts);
681:       ISDestroy(&vertices);
682:     }
683:   }
684:   if (coordSize != N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %D != %D", coordSize, N);
685:   DMLabelDestroy(&cutVertexLabel);
686:   VecRestoreArray(coordinatesLocal, &coords);
687:   VecRestoreArray(newcoords,        &ncoords);
688:   PetscObjectSetName((PetscObject) newcoords, "vertices");
689:   VecScale(newcoords, lengthScale);
690:   PetscViewerHDF5PushGroup(viewer, "/viz");
691:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
692:   PetscStackCallHDF5(H5Gclose,(groupId));
693:   PetscViewerHDF5PopGroup(viewer);
694:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
695:   VecView(newcoords, viewer);
696:   PetscViewerHDF5PopGroup(viewer);
697:   VecDestroy(&newcoords);
698:   return(0);
699: }


702: static PetscErrorCode DMPlexWriteLabels_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
703: {
704:   const PetscInt   *gpoint;
705:   PetscInt          numLabels, l;
706:   hid_t             fileId, groupId;
707:   PetscErrorCode    ierr;

710:   ISGetIndices(globalPointNumbers, &gpoint);
711:   PetscViewerHDF5PushGroup(viewer, "/labels");
712:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
713:   if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
714:   PetscViewerHDF5PopGroup(viewer);
715:   DMGetNumLabels(dm, &numLabels);
716:   for (l = 0; l < numLabels; ++l) {
717:     DMLabel         label;
718:     const char     *name;
719:     IS              valueIS, pvalueIS, globalValueIS;
720:     const PetscInt *values;
721:     PetscInt        numValues, v;
722:     PetscBool       isDepth, output;
723:     char            group[PETSC_MAX_PATH_LEN];

725:     DMGetLabelName(dm, l, &name);
726:     DMGetLabelOutput(dm, name, &output);
727:     PetscStrncmp(name, "depth", 10, &isDepth);
728:     if (isDepth || !output) continue;
729:     DMGetLabel(dm, name, &label);
730:     PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
731:     PetscViewerHDF5PushGroup(viewer, group);
732:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
733:     if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
734:     PetscViewerHDF5PopGroup(viewer);
735:     DMLabelGetValueIS(label, &valueIS);
736:     /* Must copy to a new IS on the global comm */
737:     ISGetLocalSize(valueIS, &numValues);
738:     ISGetIndices(valueIS, &values);
739:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
740:     ISRestoreIndices(valueIS, &values);
741:     ISAllGather(pvalueIS, &globalValueIS);
742:     ISDestroy(&pvalueIS);
743:     ISSortRemoveDups(globalValueIS);
744:     ISGetLocalSize(globalValueIS, &numValues);
745:     ISGetIndices(globalValueIS, &values);
746:     for (v = 0; v < numValues; ++v) {
747:       IS              stratumIS, globalStratumIS;
748:       const PetscInt *spoints = NULL;
749:       PetscInt       *gspoints, n = 0, gn, p;
750:       const char     *iname = "indices";

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

755:       if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
756:       if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
757:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
758:       PetscMalloc1(gn,&gspoints);
759:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
760:       if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
761:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
762:       if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
763:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

765:       PetscViewerHDF5PushGroup(viewer, group);
766:       ISView(globalStratumIS, viewer);
767:       PetscViewerHDF5PopGroup(viewer);
768:       ISDestroy(&globalStratumIS);
769:       ISDestroy(&stratumIS);
770:     }
771:     ISRestoreIndices(globalValueIS, &values);
772:     ISDestroy(&globalValueIS);
773:     ISDestroy(&valueIS);
774:   }
775:   ISRestoreIndices(globalPointNumbers, &gpoint);
776:   return(0);
777: }

779: /* We only write cells and vertices. Does this screw up parallel reading? */
780: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
781: {
782:   IS                globalPointNumbers;
783:   PetscViewerFormat format;
784:   PetscBool         viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
785:   PetscErrorCode    ierr;

788:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
789:   DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
790:   DMPlexWriteLabels_HDF5_Static(dm, globalPointNumbers, viewer);

792:   PetscViewerGetFormat(viewer, &format);
793:   switch (format) {
794:     case PETSC_VIEWER_HDF5_VIZ:
795:       viz_geom    = PETSC_TRUE;
796:       xdmf_topo   = PETSC_TRUE;
797:       break;
798:     case PETSC_VIEWER_HDF5_XDMF:
799:       xdmf_topo   = PETSC_TRUE;
800:       break;
801:     case PETSC_VIEWER_HDF5_PETSC:
802:       petsc_topo  = PETSC_TRUE;
803:       break;
804:     case PETSC_VIEWER_DEFAULT:
805:     case PETSC_VIEWER_NATIVE:
806:       viz_geom    = PETSC_TRUE;
807:       xdmf_topo   = PETSC_TRUE;
808:       petsc_topo  = PETSC_TRUE;
809:       break;
810:     default:
811:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
812:   }

814:   if (viz_geom)   {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
815:   if (xdmf_topo)  {DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);}
816:   if (petsc_topo) {DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);}

818:   ISDestroy(&globalPointNumbers);
819:   return(0);
820: }

822: typedef struct {
823:   PetscMPIInt rank;
824:   DM          dm;
825:   PetscViewer viewer;
826:   DMLabel     label;
827: } LabelCtx;

829: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
830: {
831:   PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
832:   DMLabel         label  = ((LabelCtx *) op_data)->label;
833:   IS              stratumIS;
834:   const PetscInt *ind;
835:   PetscInt        value, N, i;
836:   const char     *lname;
837:   char            group[PETSC_MAX_PATH_LEN];
838:   PetscErrorCode  ierr;

840:   PetscOptionsStringToInt(name, &value);
841:   ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
842:   PetscObjectSetName((PetscObject) stratumIS, "indices");
843:   PetscObjectGetName((PetscObject) label, &lname);
844:   PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
845:   PetscViewerHDF5PushGroup(viewer, group);
846:   {
847:     /* Force serial load */
848:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
849:     PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
850:     PetscLayoutSetSize(stratumIS->map, N);
851:   }
852:   ISLoad(stratumIS, viewer);
853:   PetscViewerHDF5PopGroup(viewer);
854:   ISGetLocalSize(stratumIS, &N);
855:   ISGetIndices(stratumIS, &ind);
856:   for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
857:   ISRestoreIndices(stratumIS, &ind);
858:   ISDestroy(&stratumIS);
859:   return 0;
860: }

862: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
863: {
864:   DM             dm  = ((LabelCtx *) op_data)->dm;
865:   hsize_t        idx = 0;
867:   herr_t         err;

869:   DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
870:   DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
871:   PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
872:   return err;
873: }

875: PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM dm, PetscViewer viewer)
876: {
877:   LabelCtx        ctx;
878:   hid_t           fileId, groupId;
879:   hsize_t         idx = 0;
880:   PetscErrorCode  ierr;

883:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
884:   ctx.dm     = dm;
885:   ctx.viewer = viewer;
886:   PetscViewerHDF5PushGroup(viewer, "/labels");
887:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
888:   PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
889:   PetscStackCallHDF5(H5Gclose,(groupId));
890:   PetscViewerHDF5PopGroup(viewer);
891:   return(0);
892: }

894: /* The first version will read everything onto proc 0, letting the user distribute
895:    The next will create a naive partition, and then rebalance after reading
896: */
897: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
898: {
899:   PetscSection    coordSection;
900:   Vec             coordinates;
901:   IS              orderIS, conesIS, cellsIS, orntsIS;
902:   const PetscInt *order, *cones, *cells, *ornts;
903:   PetscReal       lengthScale;
904:   PetscInt       *cone, *ornt;
905:   PetscInt        dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
906:   PetscMPIInt     rank;
907:   PetscErrorCode  ierr;

910:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
911:   /* Read toplogy */
912:   PetscViewerHDF5PushGroup(viewer, "/topology");
913:   ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
914:   PetscObjectSetName((PetscObject) orderIS, "order");
915:   ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
916:   PetscObjectSetName((PetscObject) conesIS, "cones");
917:   ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
918:   PetscObjectSetName((PetscObject) cellsIS, "cells");
919:   ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
920:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
921:   PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
922:   DMSetDimension(dm, dim);
923:   {
924:     /* Force serial load */
925:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
926:     PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
927:     PetscLayoutSetSize(orderIS->map, pEnd);
928:     PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
929:     PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
930:     PetscLayoutSetSize(conesIS->map, pEnd);
931:     pEnd = !rank ? pEnd : 0;
932:     PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
933:     PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
934:     PetscLayoutSetSize(cellsIS->map, N);
935:     PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
936:     PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
937:     PetscLayoutSetSize(orntsIS->map, N);
938:   }
939:   ISLoad(orderIS, viewer);
940:   ISLoad(conesIS, viewer);
941:   ISLoad(cellsIS, viewer);
942:   ISLoad(orntsIS, viewer);
943:   PetscViewerHDF5PopGroup(viewer);
944:   /* Read geometry */
945:   PetscViewerHDF5PushGroup(viewer, "/geometry");
946:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
947:   PetscObjectSetName((PetscObject) coordinates, "vertices");
948:   {
949:     /* Force serial load */
950:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
951:     VecSetSizes(coordinates, !rank ? N : 0, N);
952:     VecSetBlockSize(coordinates, spatialDim);
953:   }
954:   VecLoad(coordinates, viewer);
955:   PetscViewerHDF5PopGroup(viewer);
956:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
957:   VecScale(coordinates, 1.0/lengthScale);
958:   VecGetLocalSize(coordinates, &numVertices);
959:   VecGetBlockSize(coordinates, &spatialDim);
960:   numVertices /= spatialDim;
961:   /* Create Plex */
962:   DMPlexSetChart(dm, 0, pEnd);
963:   ISGetIndices(orderIS, &order);
964:   ISGetIndices(conesIS, &cones);
965:   for (p = 0; p < pEnd; ++p) {
966:     DMPlexSetConeSize(dm, order[p], cones[p]);
967:     maxConeSize = PetscMax(maxConeSize, cones[p]);
968:   }
969:   DMSetUp(dm);
970:   ISGetIndices(cellsIS, &cells);
971:   ISGetIndices(orntsIS, &ornts);
972:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
973:   for (p = 0, q = 0; p < pEnd; ++p) {
974:     for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
975:     DMPlexSetCone(dm, order[p], cone);
976:     DMPlexSetConeOrientation(dm, order[p], ornt);
977:   }
978:   PetscFree2(cone,ornt);
979:   ISRestoreIndices(orderIS, &order);
980:   ISRestoreIndices(conesIS, &cones);
981:   ISRestoreIndices(cellsIS, &cells);
982:   ISRestoreIndices(orntsIS, &ornts);
983:   ISDestroy(&orderIS);
984:   ISDestroy(&conesIS);
985:   ISDestroy(&cellsIS);
986:   ISDestroy(&orntsIS);
987:   DMPlexSymmetrize(dm);
988:   DMPlexStratify(dm);
989:   /* Create coordinates */
990:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
991:   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);
992:   DMGetCoordinateSection(dm, &coordSection);
993:   PetscSectionSetNumFields(coordSection, 1);
994:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
995:   PetscSectionSetChart(coordSection, vStart, vEnd);
996:   for (v = vStart; v < vEnd; ++v) {
997:     PetscSectionSetDof(coordSection, v, spatialDim);
998:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
999:   }
1000:   PetscSectionSetUp(coordSection);
1001:   DMSetCoordinates(dm, coordinates);
1002:   VecDestroy(&coordinates);
1003:   /* Read Labels */
1004:   DMPlexLoadLabels_HDF5_Internal(dm, viewer);
1005:   return(0);
1006: }
1007: #endif