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: 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 == 0) {
 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:   PetscViewerHDF5PushTimestepping(viewer);
 34:   PetscViewerHDF5SetTimestep(viewer, seqnum); /* seqnum < 0 jumps out above */
 35:   VecView(stamp, viewer);
 36:   PetscViewerHDF5PopTimestepping(viewer);
 37:   PetscViewerHDF5PopGroup(viewer);
 38:   VecDestroy(&stamp);
 39:   return(0);
 40: }

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

 49:   if (seqnum < 0) return(0);
 50:   MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
 51:   VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
 52:   VecSetBlockSize(stamp, 1);
 53:   PetscObjectSetName((PetscObject) stamp, seqname);
 54:   PetscViewerHDF5PushGroup(viewer, "/");
 55:   PetscViewerHDF5PushTimestepping(viewer);
 56:   PetscViewerHDF5SetTimestep(viewer, seqnum);  /* seqnum < 0 jumps out above */
 57:   VecLoad(stamp, viewer);
 58:   PetscViewerHDF5PopTimestepping(viewer);
 59:   PetscViewerHDF5PopGroup(viewer);
 60:   if (rank == 0) {
 61:     const PetscScalar *a;
 62:     PetscReal timeScale;
 63:     PetscBool istime;

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

 75: static PetscErrorCode DMPlexCreateCutVertexLabel_Private(DM dm, DMLabel cutLabel, DMLabel *cutVertexLabel)
 76: {
 77:   IS              cutcells = NULL;
 78:   const PetscInt *cutc;
 79:   PetscInt        cellHeight, vStart, vEnd, cStart, cEnd, c;
 80:   PetscErrorCode  ierr;

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

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

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

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

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

158:     PetscSectionGetNumFields(section, &numFields);
159:     DMGetLabel(dm, "periodic_cut", &cutLabel);
160:     for (f = 0; f < numFields; ++f) {
161:       Vec         subv;
162:       IS          is;
163:       const char *fname, *fgroup;
164:       char        subname[PETSC_MAX_PATH_LEN];
165:       PetscInt    pStart, pEnd;

167:       DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
168:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
169:       PetscSectionGetFieldName(section, f, &fname);
170:       if (!fname) continue;
171:       PetscViewerHDF5PushGroup(viewer, fgroup);
172:       if (cutLabel) {
173:         const PetscScalar *ga;
174:         PetscScalar       *suba;
175:         PetscInt           Nc, gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;

177:         DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
178:         PetscSectionGetFieldComponents(section, f, &Nc);
179:         for (p = pStart; p < pEnd; ++p) {
180:           PetscInt gdof, fdof = 0, val;

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

198:           PetscSectionGetDof(sectionGlobal, p, &gdof);
199:           if (gdof > 0) {
200:             PetscInt fdof, fc, f2, poff = 0;

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

245: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
246: {
247:   DM             dm;
248:   Vec            locv;
249:   PetscObject    isZero;
250:   const char    *name;
251:   PetscReal      time;

255:   VecGetDM(v, &dm);
256:   DMGetLocalVector(dm, &locv);
257:   PetscObjectGetName((PetscObject) v, &name);
258:   PetscObjectSetName((PetscObject) locv, name);
259:   PetscObjectQuery((PetscObject) v, "__Vec_bc_zero__", &isZero);
260:   PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", isZero);
261:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
262:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
263:   DMGetOutputSequenceNumber(dm, NULL, &time);
264:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
265:   PetscViewerHDF5PushGroup(viewer, "/fields");
266:   PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
267:   VecView_Plex_Local_HDF5_Internal(locv, viewer);
268:   PetscViewerPopFormat(viewer);
269:   PetscViewerHDF5PopGroup(viewer);
270:   PetscObjectCompose((PetscObject) locv, "__Vec_bc_zero__", NULL);
271:   DMRestoreLocalVector(dm, &locv);
272:   return(0);
273: }

275: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
276: {
277:   PetscBool      isseq;

281:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
282:   PetscViewerHDF5PushGroup(viewer, "/fields");
283:   if (isseq) {VecView_Seq(v, viewer);}
284:   else       {VecView_MPI(v, viewer);}
285:   PetscViewerHDF5PopGroup(viewer);
286:   return(0);
287: }

289: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
290: {
291:   DM             dm;
292:   Vec            locv;
293:   const char    *name;
294:   PetscInt       seqnum;

298:   VecGetDM(v, &dm);
299:   DMGetLocalVector(dm, &locv);
300:   PetscObjectGetName((PetscObject) v, &name);
301:   PetscObjectSetName((PetscObject) locv, name);
302:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
303:   PetscViewerHDF5PushGroup(viewer, "/fields");
304:   if (seqnum >= 0) {
305:     PetscViewerHDF5PushTimestepping(viewer);
306:     PetscViewerHDF5SetTimestep(viewer, seqnum);
307:   }
308:   VecLoad_Plex_Local(locv, viewer);
309:   if (seqnum >= 0) {
310:     PetscViewerHDF5PopTimestepping(viewer);
311:   }
312:   PetscViewerHDF5PopGroup(viewer);
313:   DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
314:   DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
315:   DMRestoreLocalVector(dm, &locv);
316:   return(0);
317: }

319: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
320: {
321:   DM             dm;
322:   PetscInt       seqnum;

326:   VecGetDM(v, &dm);
327:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
328:   PetscViewerHDF5PushGroup(viewer, "/fields");
329:   if (seqnum >= 0) {
330:     PetscViewerHDF5PushTimestepping(viewer);
331:     PetscViewerHDF5SetTimestep(viewer, seqnum);
332:   }
333:   VecLoad_Default(v, viewer);
334:   if (seqnum >= 0) {
335:     PetscViewerHDF5PopTimestepping(viewer);
336:   }
337:   PetscViewerHDF5PopGroup(viewer);
338:   return(0);
339: }

341: PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
342: {
343:   IS              orderIS, conesIS, cellsIS, orntsIS;
344:   const PetscInt *gpoint;
345:   PetscInt       *order, *sizes, *cones, *ornts;
346:   PetscInt        dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
347:   PetscErrorCode  ierr;

350:   ISGetIndices(globalPointNumbers, &gpoint);
351:   DMGetDimension(dm, &dim);
352:   DMPlexGetChart(dm, &pStart, &pEnd);
353:   for (p = pStart; p < pEnd; ++p) {
354:     if (gpoint[p] >= 0) {
355:       PetscInt coneSize;

357:       DMPlexGetConeSize(dm, p, &coneSize);
358:       conesSize += 1;
359:       cellsSize += coneSize;
360:     }
361:   }
362:   PetscMalloc1(conesSize, &order);
363:   PetscMalloc1(conesSize, &sizes);
364:   PetscMalloc1(cellsSize, &cones);
365:   PetscMalloc1(cellsSize, &ornts);
366:   for (p = pStart; p < pEnd; ++p) {
367:     if (gpoint[p] >= 0) {
368:       const PetscInt *cone, *ornt;
369:       PetscInt        coneSize, cp;

371:       DMPlexGetConeSize(dm, p, &coneSize);
372:       DMPlexGetCone(dm, p, &cone);
373:       DMPlexGetConeOrientation(dm, p, &ornt);
374:       order[s]   = gpoint[p];
375:       sizes[s++] = coneSize;
376:       for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
377:     }
378:   }
379:   if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %D != %D", s, conesSize);
380:   if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %D != %D", c, cellsSize);
381:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
382:   PetscObjectSetName((PetscObject) orderIS, "order");
383:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
384:   PetscObjectSetName((PetscObject) conesIS, "cones");
385:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
386:   PetscObjectSetName((PetscObject) cellsIS, "cells");
387:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
388:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
389:   PetscViewerHDF5PushGroup(viewer, "/topology");
390:   ISView(orderIS, viewer);
391:   ISView(conesIS, viewer);
392:   ISView(cellsIS, viewer);
393:   PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, (void *) &dim);
394:   ISView(orntsIS, viewer);
395:   PetscViewerHDF5PopGroup(viewer);
396:   ISDestroy(&orderIS);
397:   ISDestroy(&conesIS);
398:   ISDestroy(&cellsIS);
399:   ISDestroy(&orntsIS);
400:   ISRestoreIndices(globalPointNumbers, &gpoint);
401:   return(0);
402: }

404: static PetscErrorCode CreateConesIS_Private(DM dm, PetscInt cStart, PetscInt cEnd, IS globalCellNumbers, PetscInt *numCorners, IS *cellIS)
405: {
406:   PetscSF         sfPoint;
407:   DMLabel         cutLabel, cutVertexLabel = NULL;
408:   IS              globalVertexNumbers, cutvertices = NULL;
409:   const PetscInt *gcell, *gvertex, *cutverts = NULL;
410:   PetscInt       *vertices;
411:   PetscInt        conesSize = 0;
412:   PetscInt        dim, numCornersLocal = 0, cell, vStart, vEnd, vExtra = 0, v;
413:   PetscErrorCode  ierr;

416:   *numCorners = 0;
417:   DMGetDimension(dm, &dim);
418:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
419:   ISGetIndices(globalCellNumbers, &gcell);

421:   for (cell = cStart; cell < cEnd; ++cell) {
422:     PetscInt *closure = NULL;
423:     PetscInt  closureSize, v, Nc = 0;

425:     if (gcell[cell] < 0) continue;
426:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
427:     for (v = 0; v < closureSize*2; v += 2) {
428:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
429:     }
430:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
431:     conesSize += Nc;
432:     if (!numCornersLocal)           numCornersLocal = Nc;
433:     else if (numCornersLocal != Nc) numCornersLocal = 1;
434:   }
435:   MPIU_Allreduce(&numCornersLocal, numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
436:   if (numCornersLocal && (numCornersLocal != *numCorners || *numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
437:   /* Handle periodic cuts by identifying vertices which should be duplicated */
438:   DMGetLabel(dm, "periodic_cut", &cutLabel);
439:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
440:   if (cutVertexLabel) {DMLabelGetStratumIS(cutVertexLabel, 1, &cutvertices);}
441:   if (cutvertices) {
442:     ISGetIndices(cutvertices, &cutverts);
443:     ISGetLocalSize(cutvertices, &vExtra);
444:   }
445:   DMGetPointSF(dm, &sfPoint);
446:   if (cutLabel) {
447:     const PetscInt    *ilocal;
448:     const PetscSFNode *iremote;
449:     PetscInt           nroots, nleaves;

451:     PetscSFGetGraph(sfPoint, &nroots, &nleaves, &ilocal, &iremote);
452:     if (nleaves < 0) {
453:       PetscObjectReference((PetscObject) sfPoint);
454:     } else {
455:       PetscSFCreate(PetscObjectComm((PetscObject) sfPoint), &sfPoint);
456:       PetscSFSetGraph(sfPoint, nroots+vExtra, nleaves, ilocal, PETSC_USE_POINTER, iremote, PETSC_USE_POINTER);
457:     }
458:   } else {
459:     PetscObjectReference((PetscObject) sfPoint);
460:   }
461:   /* Number all vertices */
462:   DMPlexCreateNumbering_Plex(dm, vStart, vEnd+vExtra, 0, NULL, sfPoint, &globalVertexNumbers);
463:   PetscSFDestroy(&sfPoint);
464:   /* Create cones */
465:   ISGetIndices(globalVertexNumbers, &gvertex);
466:   PetscMalloc1(conesSize, &vertices);
467:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
468:     PetscInt *closure = NULL;
469:     PetscInt  closureSize, Nc = 0, p, value = -1;
470:     PetscBool replace;

472:     if (gcell[cell] < 0) continue;
473:     if (cutLabel) {DMLabelGetValue(cutLabel, cell, &value);}
474:     replace = (value == 2) ? PETSC_TRUE : PETSC_FALSE;
475:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
476:     for (p = 0; p < closureSize*2; p += 2) {
477:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
478:         closure[Nc++] = closure[p];
479:       }
480:     }
481:     DMPlexReorderCell(dm, cell, closure);
482:     for (p = 0; p < Nc; ++p) {
483:       PetscInt nv, gv = gvertex[closure[p] - vStart];

485:       if (replace) {
486:         PetscFindInt(closure[p], vExtra, cutverts, &nv);
487:         if (nv >= 0) gv = gvertex[vEnd - vStart + nv];
488:       }
489:       vertices[v++] = gv < 0 ? -(gv+1) : gv;
490:     }
491:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
492:   }
493:   ISRestoreIndices(globalVertexNumbers, &gvertex);
494:   ISDestroy(&globalVertexNumbers);
495:   ISRestoreIndices(globalCellNumbers, &gcell);
496:   if (cutvertices) {ISRestoreIndices(cutvertices, &cutverts);}
497:   ISDestroy(&cutvertices);
498:   DMLabelDestroy(&cutVertexLabel);
499:   if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %D != %D", v, conesSize);
500:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, cellIS);
501:   PetscLayoutSetBlockSize((*cellIS)->map, *numCorners);
502:   PetscObjectSetName((PetscObject) *cellIS, "cells");
503:   return(0);
504: }

506: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, IS globalCellNumbers, PetscViewer viewer)
507: {
508:   DM              cdm;
509:   DMLabel         depthLabel, ctLabel;
510:   IS              cellIS;
511:   PetscInt        dim, depth, cellHeight, c;
512:   hid_t           fileId, groupId;
513:   PetscErrorCode  ierr;

516:   PetscViewerHDF5PushGroup(viewer, "/viz");
517:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
518:   PetscStackCallHDF5(H5Gclose,(groupId));

520:   PetscViewerHDF5PopGroup(viewer);
521:   DMGetDimension(dm, &dim);
522:   DMPlexGetDepth(dm, &depth);
523:   DMGetCoordinateDM(dm, &cdm);
524:   DMPlexGetVTKCellHeight(dm, &cellHeight);
525:   DMPlexGetDepthLabel(dm, &depthLabel);
526:   DMPlexGetCellTypeLabel(dm, &ctLabel);
527:   for (c = 0; c < DM_NUM_POLYTOPES; ++c) {
528:     const DMPolytopeType ict = (DMPolytopeType) c;
529:     PetscInt             pStart, pEnd, dep, numCorners, n = 0;
530:     PetscBool            output = PETSC_FALSE, doOutput;

532:     if (ict == DM_POLYTOPE_FV_GHOST) continue;
533:     DMLabelGetStratumBounds(ctLabel, ict, &pStart, &pEnd);
534:     if (pStart >= 0) {
535:       DMLabelGetValue(depthLabel, pStart, &dep);
536:       if (dep == depth - cellHeight) output = PETSC_TRUE;
537:     }
538:     MPI_Allreduce(&output, &doOutput, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
539:     if (!doOutput) continue;
540:     CreateConesIS_Private(dm, pStart, pEnd, globalCellNumbers, &numCorners,  &cellIS);
541:     if (!n) {
542:       PetscViewerHDF5PushGroup(viewer, "/viz/topology");
543:     } else {
544:       char group[PETSC_MAX_PATH_LEN];

546:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/viz/topology_%D", n);
547:       PetscViewerHDF5PushGroup(viewer, group);
548:     }
549:     ISView(cellIS, viewer);
550:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_corners", PETSC_INT, (void *) &numCorners);
551:     PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) cellIS, "cell_dim",     PETSC_INT, (void *) &dim);
552:     ISDestroy(&cellIS);
553:     PetscViewerHDF5PopGroup(viewer);
554:     ++n;
555:   }
556:   return(0);
557: }

559: PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM dm, PetscViewer viewer)
560: {
561:   DM             cdm;
562:   Vec            coordinates, newcoords;
563:   PetscReal      lengthScale;
564:   PetscInt       m, M, bs;

568:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
569:   DMGetCoordinateDM(dm, &cdm);
570:   DMGetCoordinates(dm, &coordinates);
571:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
572:   PetscObjectSetName((PetscObject) newcoords, "vertices");
573:   VecGetSize(coordinates, &M);
574:   VecGetLocalSize(coordinates, &m);
575:   VecSetSizes(newcoords, m, M);
576:   VecGetBlockSize(coordinates, &bs);
577:   VecSetBlockSize(newcoords, bs);
578:   VecSetType(newcoords,VECSTANDARD);
579:   VecCopy(coordinates, newcoords);
580:   VecScale(newcoords, lengthScale);
581:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
582:   PetscViewerHDF5PushGroup(viewer, "/geometry");
583:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
584:   VecView(newcoords, viewer);
585:   PetscViewerPopFormat(viewer);
586:   PetscViewerHDF5PopGroup(viewer);
587:   VecDestroy(&newcoords);
588:   return(0);
589: }

591: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
592: {
593:   DM               cdm;
594:   Vec              coordinatesLocal, newcoords;
595:   PetscSection     cSection, cGlobalSection;
596:   PetscScalar     *coords, *ncoords;
597:   DMLabel          cutLabel, cutVertexLabel = NULL;
598:   const PetscReal *L;
599:   const DMBoundaryType *bd;
600:   PetscReal        lengthScale;
601:   PetscInt         vStart, vEnd, v, bs, N, coordSize, dof, off, d;
602:   PetscBool        localized, embedded;
603:   hid_t            fileId, groupId;
604:   PetscErrorCode   ierr;

607:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
608:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
609:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
610:   VecGetBlockSize(coordinatesLocal, &bs);
611:   DMGetCoordinatesLocalized(dm, &localized);
612:   if (localized == PETSC_FALSE) return(0);
613:   DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
614:   DMGetCoordinateDM(dm, &cdm);
615:   DMGetLocalSection(cdm, &cSection);
616:   DMGetGlobalSection(cdm, &cGlobalSection);
617:   DMGetLabel(dm, "periodic_cut", &cutLabel);
618:   N    = 0;

620:   DMPlexCreateCutVertexLabel_Private(dm, cutLabel, &cutVertexLabel);
621:   VecCreate(PetscObjectComm((PetscObject) dm), &newcoords);
622:   PetscSectionGetDof(cSection, vStart, &dof);
623:   PetscPrintf(PETSC_COMM_SELF, "DOF: %D\n", dof);
624:   embedded  = (PetscBool) (L && dof == 2 && !cutLabel);
625:   if (cutVertexLabel) {
626:     DMLabelGetStratumSize(cutVertexLabel, 1, &v);
627:     N   += dof*v;
628:   }
629:   for (v = vStart; v < vEnd; ++v) {
630:     PetscSectionGetDof(cGlobalSection, v, &dof);
631:     if (dof < 0) continue;
632:     if (embedded) N += dof+1;
633:     else          N += dof;
634:   }
635:   if (embedded) {VecSetBlockSize(newcoords, bs+1);}
636:   else          {VecSetBlockSize(newcoords, bs);}
637:   VecSetSizes(newcoords, N, PETSC_DETERMINE);
638:   VecSetType(newcoords, VECSTANDARD);
639:   VecGetArray(coordinatesLocal, &coords);
640:   VecGetArray(newcoords,        &ncoords);
641:   coordSize = 0;
642:   for (v = vStart; v < vEnd; ++v) {
643:     PetscSectionGetDof(cGlobalSection, v, &dof);
644:     PetscSectionGetOffset(cSection, v, &off);
645:     if (dof < 0) continue;
646:     if (embedded) {
647:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
648:         PetscReal theta, phi, r, R;
649:         /* XY-periodic */
650:         /* Suppose its an y-z circle, then
651:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
652:            and the circle in that plane is
653:              \hat r cos(phi) + \hat x sin(phi) */
654:         theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
655:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
656:         r     = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
657:         R     = L[1]/(2.0*PETSC_PI);
658:         ncoords[coordSize++] =  PetscSinReal(phi) * r;
659:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
660:         ncoords[coordSize++] =  PetscSinReal(theta) * (R + r * PetscCosReal(phi));
661:       } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
662:         /* X-periodic */
663:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
664:         ncoords[coordSize++] = coords[off+1];
665:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
666:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
667:         /* Y-periodic */
668:         ncoords[coordSize++] = coords[off+0];
669:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
670:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
671:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
672:         PetscReal phi, r, R;
673:         /* Mobius strip */
674:         /* Suppose its an x-z circle, then
675:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
676:            and in that plane we rotate by pi as we go around the circle
677:              \hat r cos(phi/2) + \hat y sin(phi/2) */
678:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
679:         R     = L[0];
680:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
681:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
682:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
683:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
684:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
685:     } else {
686:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
687:     }
688:   }
689:   if (cutVertexLabel) {
690:     IS              vertices;
691:     const PetscInt *verts;
692:     PetscInt        n;

694:     DMLabelGetStratumIS(cutVertexLabel, 1, &vertices);
695:     if (vertices) {
696:       ISGetIndices(vertices, &verts);
697:       ISGetLocalSize(vertices, &n);
698:       for (v = 0; v < n; ++v) {
699:         PetscSectionGetDof(cSection, verts[v], &dof);
700:         PetscSectionGetOffset(cSection, verts[v], &off);
701:         for (d = 0; d < dof; ++d) ncoords[coordSize++] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
702:       }
703:       ISRestoreIndices(vertices, &verts);
704:       ISDestroy(&vertices);
705:     }
706:   }
707:   if (coordSize != N) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched sizes: %D != %D", coordSize, N);
708:   DMLabelDestroy(&cutVertexLabel);
709:   VecRestoreArray(coordinatesLocal, &coords);
710:   VecRestoreArray(newcoords,        &ncoords);
711:   PetscObjectSetName((PetscObject) newcoords, "vertices");
712:   VecScale(newcoords, lengthScale);
713:   PetscViewerHDF5PushGroup(viewer, "/viz");
714:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
715:   PetscStackCallHDF5(H5Gclose,(groupId));
716:   PetscViewerHDF5PopGroup(viewer);
717:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
718:   VecView(newcoords, viewer);
719:   PetscViewerHDF5PopGroup(viewer);
720:   VecDestroy(&newcoords);
721:   return(0);
722: }

724: PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM dm, IS globalPointNumbers, PetscViewer viewer)
725: {
726:   const PetscInt   *gpoint;
727:   PetscInt          numLabels, l;
728:   hid_t             fileId, groupId;
729:   PetscErrorCode    ierr;

732:   ISGetIndices(globalPointNumbers, &gpoint);
733:   PetscViewerHDF5PushGroup(viewer, "/labels");
734:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
735:   if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
736:   PetscViewerHDF5PopGroup(viewer);
737:   DMGetNumLabels(dm, &numLabels);
738:   for (l = 0; l < numLabels; ++l) {
739:     DMLabel         label;
740:     const char     *name;
741:     IS              valueIS, pvalueIS, globalValueIS;
742:     const PetscInt *values;
743:     PetscInt        numValues, v;
744:     PetscBool       isDepth, output;
745:     char            group[PETSC_MAX_PATH_LEN];

747:     DMGetLabelName(dm, l, &name);
748:     DMGetLabelOutput(dm, name, &output);
749:     PetscStrncmp(name, "depth", 10, &isDepth);
750:     if (isDepth || !output) continue;
751:     DMGetLabel(dm, name, &label);
752:     PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
753:     PetscViewerHDF5PushGroup(viewer, group);
754:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
755:     if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
756:     PetscViewerHDF5PopGroup(viewer);
757:     DMLabelGetValueIS(label, &valueIS);
758:     /* Must copy to a new IS on the global comm */
759:     ISGetLocalSize(valueIS, &numValues);
760:     ISGetIndices(valueIS, &values);
761:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
762:     ISRestoreIndices(valueIS, &values);
763:     ISAllGather(pvalueIS, &globalValueIS);
764:     ISDestroy(&pvalueIS);
765:     ISSortRemoveDups(globalValueIS);
766:     ISGetLocalSize(globalValueIS, &numValues);
767:     ISGetIndices(globalValueIS, &values);
768:     for (v = 0; v < numValues; ++v) {
769:       IS              stratumIS, globalStratumIS;
770:       const PetscInt *spoints = NULL;
771:       PetscInt       *gspoints, n = 0, gn, p;
772:       const char     *iname = "indices";

774:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%D", name, values[v]);
775:       DMLabelGetStratumIS(label, values[v], &stratumIS);

777:       if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
778:       if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
779:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
780:       PetscMalloc1(gn,&gspoints);
781:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
782:       if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
783:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
784:       if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
785:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

787:       PetscViewerHDF5PushGroup(viewer, group);
788:       ISView(globalStratumIS, viewer);
789:       PetscViewerHDF5PopGroup(viewer);
790:       ISDestroy(&globalStratumIS);
791:       ISDestroy(&stratumIS);
792:     }
793:     ISRestoreIndices(globalValueIS, &values);
794:     ISDestroy(&globalValueIS);
795:     ISDestroy(&valueIS);
796:   }
797:   ISRestoreIndices(globalPointNumbers, &gpoint);
798:   return(0);
799: }

801: /* We only write cells and vertices. Does this screw up parallel reading? */
802: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
803: {
804:   IS                globalPointNumbers;
805:   PetscViewerFormat format;
806:   PetscBool         viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
807:   PetscErrorCode    ierr;

810:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
811:   DMPlexCoordinatesView_HDF5_Internal(dm, viewer);
812:   DMPlexLabelsView_HDF5_Internal(dm, globalPointNumbers, viewer);

814:   PetscViewerHDF5WriteAttribute(viewer, NULL, "petsc_version_git", PETSC_STRING, PETSC_VERSION_GIT);

816:   PetscViewerGetFormat(viewer, &format);
817:   switch (format) {
818:     case PETSC_VIEWER_HDF5_VIZ:
819:       viz_geom    = PETSC_TRUE;
820:       xdmf_topo   = PETSC_TRUE;
821:       break;
822:     case PETSC_VIEWER_HDF5_XDMF:
823:       xdmf_topo   = PETSC_TRUE;
824:       break;
825:     case PETSC_VIEWER_HDF5_PETSC:
826:       petsc_topo  = PETSC_TRUE;
827:       break;
828:     case PETSC_VIEWER_DEFAULT:
829:     case PETSC_VIEWER_NATIVE:
830:       viz_geom    = PETSC_TRUE;
831:       xdmf_topo   = PETSC_TRUE;
832:       petsc_topo  = PETSC_TRUE;
833:       break;
834:     default:
835:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
836:   }

838:   if (viz_geom)   {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
839:   if (xdmf_topo)  {DMPlexWriteTopology_Vertices_HDF5_Static(dm, globalPointNumbers, viewer);}
840:   if (petsc_topo) {DMPlexTopologyView_HDF5_Internal(dm, globalPointNumbers, viewer);}

842:   ISDestroy(&globalPointNumbers);
843:   return(0);
844: }

846: PetscErrorCode DMPlexSectionView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm)
847: {
848:   MPI_Comm       comm;
849:   const char    *topologydm_name;
850:   const char    *sectiondm_name;
851:   PetscSection   gsection;

855:   PetscObjectGetComm((PetscObject)sectiondm, &comm);
856:   PetscObjectGetName((PetscObject)dm, &topologydm_name);
857:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
858:   PetscViewerHDF5PushGroup(viewer, "topologies");
859:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
860:   PetscViewerHDF5PushGroup(viewer, "dms");
861:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
862:   DMGetGlobalSection(sectiondm, &gsection);
863:   /* Save raw section */
864:   PetscSectionView(gsection, viewer);
865:   /* Save plex wrapper */
866:   {
867:     PetscInt        pStart, pEnd, p, n;
868:     IS              globalPointNumbers;
869:     const PetscInt *gpoints;
870:     IS              orderIS;
871:     PetscInt       *order;

873:     PetscSectionGetChart(gsection, &pStart, &pEnd);
874:     DMPlexCreatePointNumbering(dm, &globalPointNumbers);
875:     ISGetIndices(globalPointNumbers, &gpoints);
876:     for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) n++;
877:     /* "order" is an array of global point numbers.
878:        When loading, it is used with topology/order array
879:        to match section points with plex topology points. */
880:     PetscMalloc1(n, &order);
881:     for (p = pStart, n = 0; p < pEnd; ++p) if (gpoints[p] >= 0) order[n++] = gpoints[p];
882:     ISRestoreIndices(globalPointNumbers, &gpoints);
883:     ISDestroy(&globalPointNumbers);
884:     ISCreateGeneral(comm, n, order, PETSC_OWN_POINTER, &orderIS);
885:     PetscObjectSetName((PetscObject)orderIS, "order");
886:     ISView(orderIS, viewer);
887:     ISDestroy(&orderIS);
888:   }
889:   PetscViewerHDF5PopGroup(viewer);
890:   PetscViewerHDF5PopGroup(viewer);
891:   PetscViewerHDF5PopGroup(viewer);
892:   PetscViewerHDF5PopGroup(viewer);
893:   return(0);
894: }

896: PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
897: {
898:   const char     *topologydm_name;
899:   const char     *sectiondm_name;
900:   const char     *vec_name;
901:   PetscInt        bs;
902:   PetscErrorCode  ierr;

905:   /* Check consistency */
906:   {
907:     PetscSF   pointsf, pointsf1;

909:     DMGetPointSF(dm, &pointsf);
910:     DMGetPointSF(sectiondm, &pointsf1);
911:     if (pointsf1 != pointsf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
912:   }
913:   PetscObjectGetName((PetscObject)dm, &topologydm_name);
914:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
915:   PetscObjectGetName((PetscObject)vec, &vec_name);
916:   PetscViewerHDF5PushGroup(viewer, "topologies");
917:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
918:   PetscViewerHDF5PushGroup(viewer, "dms");
919:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
920:   PetscViewerHDF5PushGroup(viewer, "vecs");
921:   PetscViewerHDF5PushGroup(viewer, vec_name);
922:   VecGetBlockSize(vec, &bs);
923:   PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
924:   VecSetBlockSize(vec, 1);
925:   /* VecView(vec, viewer) would call (*vec->opt->view)(vec, viewer), but,    */
926:   /* if vec was created with DMGet{Global, Local}Vector(), vec->opt->view    */
927:   /* is set to VecView_Plex, which would save vec in a predefined location.  */
928:   /* To save vec in where we want, we create a new Vec (temp) with           */
929:   /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
930:   {
931:     Vec                temp;
932:     const PetscScalar *array;
933:     PetscLayout        map;

935:     VecCreate(PetscObjectComm((PetscObject)vec), &temp);
936:     PetscObjectSetName((PetscObject)temp, vec_name);
937:     VecGetLayout(vec, &map);
938:     VecSetLayout(temp, map);
939:     VecSetUp(temp);
940:     VecGetArrayRead(vec, &array);
941:     VecPlaceArray(temp, array);
942:     VecView(temp, viewer);
943:     VecResetArray(temp);
944:     VecRestoreArrayRead(vec, &array);
945:     VecDestroy(&temp);
946:   }
947:   VecSetBlockSize(vec, bs);
948:   PetscViewerHDF5PopGroup(viewer);
949:   PetscViewerHDF5PopGroup(viewer);
950:   PetscViewerHDF5PopGroup(viewer);
951:   PetscViewerHDF5PopGroup(viewer);
952:   PetscViewerHDF5PopGroup(viewer);
953:   PetscViewerHDF5PopGroup(viewer);
954:   return(0);
955: }

957: PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, Vec vec)
958: {
959:   MPI_Comm        comm;
960:   const char     *topologydm_name;
961:   const char     *sectiondm_name;
962:   const char     *vec_name;
963:   PetscSection    section;
964:   PetscBool       includesConstraints;
965:   Vec             gvec;
966:   PetscInt        m, bs;
967:   PetscErrorCode  ierr;

970:   PetscObjectGetComm((PetscObject)dm, &comm);
971:   /* Check consistency */
972:   {
973:     PetscSF   pointsf, pointsf1;

975:     DMGetPointSF(dm, &pointsf);
976:     DMGetPointSF(sectiondm, &pointsf1);
977:     if (pointsf1 != pointsf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
978:   }
979:   PetscObjectGetName((PetscObject)dm, &topologydm_name);
980:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
981:   PetscObjectGetName((PetscObject)vec, &vec_name);
982:   PetscViewerHDF5PushGroup(viewer, "topologies");
983:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
984:   PetscViewerHDF5PushGroup(viewer, "dms");
985:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
986:   PetscViewerHDF5PushGroup(viewer, "vecs");
987:   PetscViewerHDF5PushGroup(viewer, vec_name);
988:   VecGetBlockSize(vec, &bs);
989:   PetscViewerHDF5WriteAttribute(viewer, NULL, "blockSize", PETSC_INT, (void *) &bs);
990:   VecCreate(comm, &gvec);
991:   PetscObjectSetName((PetscObject)gvec, vec_name);
992:   DMGetGlobalSection(sectiondm, &section);
993:   PetscSectionGetIncludesConstraints(section, &includesConstraints);
994:   if (includesConstraints) {PetscSectionGetStorageSize(section, &m);}
995:   else {PetscSectionGetConstrainedStorageSize(section, &m);}
996:   VecSetSizes(gvec, m, PETSC_DECIDE);
997:   VecSetUp(gvec);
998:   DMLocalToGlobalBegin(sectiondm, vec, INSERT_VALUES, gvec);
999:   DMLocalToGlobalEnd(sectiondm, vec, INSERT_VALUES, gvec);
1000:   VecView(gvec, viewer);
1001:   VecDestroy(&gvec);
1002:   PetscViewerHDF5PopGroup(viewer);
1003:   PetscViewerHDF5PopGroup(viewer);
1004:   PetscViewerHDF5PopGroup(viewer);
1005:   PetscViewerHDF5PopGroup(viewer);
1006:   PetscViewerHDF5PopGroup(viewer);
1007:   PetscViewerHDF5PopGroup(viewer);
1008:   return(0);
1009: }

1011: typedef struct {
1012:   PetscMPIInt rank;
1013:   DM          dm;
1014:   PetscViewer viewer;
1015:   DMLabel     label;
1016: } LabelCtx;

1018: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
1019: {
1020:   PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
1021:   DMLabel         label  = ((LabelCtx *) op_data)->label;
1022:   IS              stratumIS;
1023:   const PetscInt *ind;
1024:   PetscInt        value, N, i;
1025:   const char     *lname;
1026:   char            group[PETSC_MAX_PATH_LEN];
1027:   PetscErrorCode  ierr;

1029:   PetscOptionsStringToInt(name, &value);
1030:   ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
1031:   PetscObjectSetName((PetscObject) stratumIS, "indices");
1032:   PetscObjectGetName((PetscObject) label, &lname);
1033:   PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
1034:   PetscViewerHDF5PushGroup(viewer, group);
1035:   {
1036:     /* Force serial load */
1037:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
1038:     PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
1039:     PetscLayoutSetSize(stratumIS->map, N);
1040:   }
1041:   ISLoad(stratumIS, viewer);
1042:   PetscViewerHDF5PopGroup(viewer);
1043:   ISGetLocalSize(stratumIS, &N);
1044:   ISGetIndices(stratumIS, &ind);
1045:   for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
1046:   ISRestoreIndices(stratumIS, &ind);
1047:   ISDestroy(&stratumIS);
1048:   return 0;
1049: }

1051: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
1052: {
1053:   DM             dm  = ((LabelCtx *) op_data)->dm;
1054:   hsize_t        idx = 0;
1056:   herr_t         err;

1058:   DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
1059:   DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
1060:   PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
1061:   return err;
1062: }

1064: PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1065: {
1066:   LabelCtx        ctx;
1067:   hid_t           fileId, groupId;
1068:   hsize_t         idx = 0;
1069:   PetscErrorCode  ierr;

1072:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
1073:   ctx.dm     = dm;
1074:   ctx.viewer = viewer;
1075:   PetscViewerHDF5PushGroup(viewer, "/labels");
1076:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
1077:   PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
1078:   PetscStackCallHDF5(H5Gclose,(groupId));
1079:   PetscViewerHDF5PopGroup(viewer);
1080:   return(0);
1081: }

1083: PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM dm, PetscViewer viewer, PetscSF *sf)
1084: {
1085:   MPI_Comm        comm;
1086:   IS              orderIS, conesIS, cellsIS, orntsIS;
1087:   const PetscInt *order, *cones, *cells, *ornts;
1088:   PetscInt       *cone, *ornt;
1089:   PetscInt        dim, N, Np, pEnd, p, q, maxConeSize = 0, c;
1090:   PetscMPIInt     size, rank;
1091:   PetscErrorCode  ierr;

1094:   PetscObjectGetComm((PetscObject)dm, &comm);
1095:   MPI_Comm_size(comm, &size);
1096:   MPI_Comm_rank(comm, &rank);
1097:   /* Read toplogy */
1098:   PetscViewerHDF5PushGroup(viewer, "/topology");
1099:   ISCreate(comm, &orderIS);
1100:   PetscObjectSetName((PetscObject) orderIS, "order");
1101:   ISCreate(comm, &conesIS);
1102:   PetscObjectSetName((PetscObject) conesIS, "cones");
1103:   ISCreate(comm, &cellsIS);
1104:   PetscObjectSetName((PetscObject) cellsIS, "cells");
1105:   ISCreate(comm, &orntsIS);
1106:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
1107:   PetscViewerHDF5ReadObjectAttribute(viewer, (PetscObject) cellsIS, "cell_dim", PETSC_INT, NULL, &dim);
1108:   DMSetDimension(dm, dim);
1109:   {
1110:     /* Force serial load */
1111:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &Np);
1112:     PetscLayoutSetLocalSize(orderIS->map, rank == 0 ? Np : 0);
1113:     PetscLayoutSetSize(orderIS->map, Np);
1114:     PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &Np);
1115:     PetscLayoutSetLocalSize(conesIS->map, rank == 0 ? Np : 0);
1116:     PetscLayoutSetSize(conesIS->map, Np);
1117:     pEnd = rank == 0 ? Np : 0;
1118:     PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
1119:     PetscLayoutSetLocalSize(cellsIS->map, rank == 0 ? N : 0);
1120:     PetscLayoutSetSize(cellsIS->map, N);
1121:     PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
1122:     PetscLayoutSetLocalSize(orntsIS->map, rank == 0 ? N : 0);
1123:     PetscLayoutSetSize(orntsIS->map, N);
1124:   }
1125:   ISLoad(orderIS, viewer);
1126:   ISLoad(conesIS, viewer);
1127:   ISLoad(cellsIS, viewer);
1128:   ISLoad(orntsIS, viewer);
1129:   PetscViewerHDF5PopGroup(viewer);
1130:   /* Create Plex */
1131:   DMPlexSetChart(dm, 0, pEnd);
1132:   ISGetIndices(orderIS, &order);
1133:   ISGetIndices(conesIS, &cones);
1134:   for (p = 0; p < pEnd; ++p) {
1135:     DMPlexSetConeSize(dm, order[p], cones[p]);
1136:     maxConeSize = PetscMax(maxConeSize, cones[p]);
1137:   }
1138:   DMSetUp(dm);
1139:   ISGetIndices(cellsIS, &cells);
1140:   ISGetIndices(orntsIS, &ornts);
1141:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
1142:   for (p = 0, q = 0; p < pEnd; ++p) {
1143:     for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
1144:     DMPlexSetCone(dm, order[p], cone);
1145:     DMPlexSetConeOrientation(dm, order[p], ornt);
1146:   }
1147:   PetscFree2(cone,ornt);
1148:   /* Create global section migration SF */
1149:   if (sf) {
1150:     PetscLayout  layout;
1151:     PetscInt    *globalIndices;

1153:     PetscMalloc1(pEnd, &globalIndices);
1154:     /* plex point == globalPointNumber in this case */
1155:     for (p = 0; p < pEnd; ++p) globalIndices[p] = p;
1156:     PetscLayoutCreate(comm, &layout);
1157:     PetscLayoutSetSize(layout, Np);
1158:     PetscLayoutSetBlockSize(layout, 1);
1159:     PetscLayoutSetUp(layout);
1160:     PetscSFCreate(comm, sf);
1161:     PetscSFSetFromOptions(*sf);
1162:     PetscSFSetGraphLayout(*sf, layout, pEnd, NULL, PETSC_OWN_POINTER, globalIndices);
1163:     PetscLayoutDestroy(&layout);
1164:     PetscFree(globalIndices);
1165:   }
1166:   /*  */
1167:   ISRestoreIndices(orderIS, &order);
1168:   ISRestoreIndices(conesIS, &cones);
1169:   ISRestoreIndices(cellsIS, &cells);
1170:   ISRestoreIndices(orntsIS, &ornts);
1171:   ISDestroy(&orderIS);
1172:   ISDestroy(&conesIS);
1173:   ISDestroy(&cellsIS);
1174:   ISDestroy(&orntsIS);
1175:   DMPlexSymmetrize(dm);
1176:   DMPlexStratify(dm);
1177:   return(0);
1178: }

1180: PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1181: {
1182:   PetscSection    coordSection;
1183:   Vec             coordinates;
1184:   PetscReal       lengthScale;
1185:   PetscInt        spatialDim, N, numVertices, vStart, vEnd, v;
1186:   PetscMPIInt     rank;
1187:   PetscErrorCode  ierr;

1190:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
1191:   /* Read geometry */
1192:   PetscViewerHDF5PushGroup(viewer, "/geometry");
1193:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
1194:   PetscObjectSetName((PetscObject) coordinates, "vertices");
1195:   {
1196:     /* Force serial load */
1197:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
1198:     VecSetSizes(coordinates, !rank ? N : 0, N);
1199:     VecSetBlockSize(coordinates, spatialDim);
1200:   }
1201:   VecLoad(coordinates, viewer);
1202:   PetscViewerHDF5PopGroup(viewer);
1203:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
1204:   VecScale(coordinates, 1.0/lengthScale);
1205:   VecGetLocalSize(coordinates, &numVertices);
1206:   VecGetBlockSize(coordinates, &spatialDim);
1207:   numVertices /= spatialDim;
1208:   /* Create coordinates */
1209:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1210:   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);
1211:   DMGetCoordinateSection(dm, &coordSection);
1212:   PetscSectionSetNumFields(coordSection, 1);
1213:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
1214:   PetscSectionSetChart(coordSection, vStart, vEnd);
1215:   for (v = vStart; v < vEnd; ++v) {
1216:     PetscSectionSetDof(coordSection, v, spatialDim);
1217:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
1218:   }
1219:   PetscSectionSetUp(coordSection);
1220:   DMSetCoordinates(dm, coordinates);
1221:   VecDestroy(&coordinates);
1222:   return(0);
1223: }

1225: /* The first version will read everything onto proc 0, letting the user distribute
1226:    The next will create a naive partition, and then rebalance after reading
1227: */
1228: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
1229: {
1230:   PetscErrorCode  ierr;

1233:   DMPlexTopologyLoad_HDF5_Internal(dm, viewer, NULL);
1234:   DMPlexCoordinatesLoad_HDF5_Internal(dm, viewer);
1235:   DMPlexLabelsLoad_HDF5_Internal(dm, viewer);
1236:   return(0);
1237: }

1239: static PetscErrorCode DMPlexSectionLoad_HDF5_Internal_CreateDataSF(PetscSection rootSection, PetscLayout layout, PetscInt globalOffsets[], PetscSection leafSection, PetscSF *sectionSF)
1240: {
1241:   MPI_Comm        comm;
1242:   PetscInt        pStart, pEnd, p, m;
1243:   PetscInt       *goffs, *ilocal;
1244:   PetscBool       rootIncludeConstraints, leafIncludeConstraints;
1245:   PetscErrorCode  ierr;

1248:   PetscObjectGetComm((PetscObject)leafSection, &comm);
1249:   PetscSectionGetChart(leafSection, &pStart, &pEnd);
1250:   PetscSectionGetIncludesConstraints(rootSection, &rootIncludeConstraints);
1251:   PetscSectionGetIncludesConstraints(leafSection, &leafIncludeConstraints);
1252:   if (rootIncludeConstraints && leafIncludeConstraints) {PetscSectionGetStorageSize(leafSection, &m);}
1253:   else {PetscSectionGetConstrainedStorageSize(leafSection, &m);}
1254:   PetscMalloc1(m, &ilocal);
1255:   PetscMalloc1(m, &goffs);
1256:   /* Currently, PetscSFDistributeSection() returns globalOffsets[] only */
1257:   /* for the top-level section (not for each field), so one must have   */
1258:   /* rootSection->pointMajor == PETSC_TRUE.                             */
1259:   if (!rootSection->pointMajor) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for field major ordering");
1260:   /* Currently, we also assume that leafSection->pointMajor == PETSC_TRUE. */
1261:   if (!leafSection->pointMajor) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for field major ordering");
1262:   for (p = pStart, m = 0; p < pEnd; ++p) {
1263:     PetscInt        dof, cdof, i, j, off, goff;
1264:     const PetscInt *cinds;

1266:     PetscSectionGetDof(leafSection, p, &dof);
1267:     if (dof < 0) continue;
1268:     goff = globalOffsets[p-pStart];
1269:     PetscSectionGetOffset(leafSection, p, &off);
1270:     PetscSectionGetConstraintDof(leafSection, p, &cdof);
1271:     PetscSectionGetConstraintIndices(leafSection, p, &cinds);
1272:     for (i = 0, j = 0; i < dof; ++i) {
1273:       PetscBool constrained = (PetscBool) (j < cdof && i == cinds[j]);

1275:       if (!constrained || (leafIncludeConstraints && rootIncludeConstraints)) {ilocal[m] = off++; goffs[m++] = goff++;}
1276:       else if (leafIncludeConstraints && !rootIncludeConstraints) ++off;
1277:       else if (!leafIncludeConstraints &&  rootIncludeConstraints) ++goff;
1278:       if (constrained) ++j;
1279:     }
1280:   }
1281:   PetscSFCreate(comm, sectionSF);
1282:   PetscSFSetFromOptions(*sectionSF);
1283:   PetscSFSetGraphLayout(*sectionSF, layout, m, ilocal, PETSC_OWN_POINTER, goffs);
1284:   PetscFree(goffs);
1285:   return(0);
1286: }

1288: PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sfXB, PetscSF *gsf, PetscSF *lsf)
1289: {
1290:   MPI_Comm       comm;
1291:   PetscMPIInt    size, rank;
1292:   const char    *topologydm_name;
1293:   const char    *sectiondm_name;
1294:   PetscSection   sectionA, sectionB;
1295:   PetscInt       NX, nX, n, i;
1296:   PetscSF        sfAB;

1300:   PetscObjectGetComm((PetscObject)dm, &comm);
1301:   MPI_Comm_size(comm, &size);
1302:   MPI_Comm_rank(comm, &rank);
1303:   PetscObjectGetName((PetscObject)dm, &topologydm_name);
1304:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1305:   PetscViewerHDF5PushGroup(viewer, "topologies");
1306:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1307:   PetscViewerHDF5ReadSizes(viewer, "/topology/order", NULL, &NX);
1308:   PetscViewerHDF5PushGroup(viewer, "dms");
1309:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1310:   /* A: on-disk points                        */
1311:   /* X: list of global point numbers, [0, NX) */
1312:   /* B: plex points                           */
1313:   /* Load raw section (sectionA)              */
1314:   PetscSectionCreate(comm, &sectionA);
1315:   PetscSectionLoad(sectionA, viewer);
1316:   PetscSectionGetChart(sectionA, NULL, &n);
1317:   /* Create sfAB: A -> B */
1318: #if defined(PETSC_USE_DEBUG)
1319:   {
1320:     PetscInt  N, N1;

1322:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &N1);
1323:     MPI_Allreduce(&n, &N, 1, MPIU_INT, MPI_SUM, comm);
1324:     if (N1 != N) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Mismatching sizes: on-disk order array size (%D) != number of loaded section points (%D)", N1, N);
1325:   }
1326: #endif
1327:   {
1328:     IS              orderIS;
1329:     const PetscInt *gpoints;
1330:     PetscSF         sfXA, sfAX;
1331:     PetscLayout     layout;
1332:     PetscSFNode    *owners, *buffer;
1333:     PetscInt        nleaves;
1334:     PetscInt       *ilocal;
1335:     PetscSFNode    *iremote;

1337:     /* Create sfAX: A -> X */
1338:     ISCreate(comm, &orderIS);
1339:     PetscObjectSetName((PetscObject)orderIS, "order");
1340:     PetscLayoutSetLocalSize(orderIS->map, n);
1341:     ISLoad(orderIS, viewer);
1342:     PetscLayoutCreate(comm, &layout);
1343:     PetscLayoutSetSize(layout, NX);
1344:     PetscLayoutSetBlockSize(layout, 1);
1345:     PetscLayoutSetUp(layout);
1346:     PetscSFCreate(comm, &sfXA);
1347:     ISGetIndices(orderIS, &gpoints);
1348:     PetscSFSetGraphLayout(sfXA, layout, n, NULL, PETSC_OWN_POINTER, gpoints);
1349:     ISRestoreIndices(orderIS, &gpoints);
1350:     ISDestroy(&orderIS);
1351:     PetscLayoutDestroy(&layout);
1352:     PetscSFGetGraph(sfXA, &nX, NULL, NULL, NULL);
1353:     PetscMalloc1(n, &owners);
1354:     PetscMalloc1(nX, &buffer);
1355:     for (i = 0; i < n; ++i) {owners[i].rank = rank; owners[i].index = i;}
1356:     for (i = 0; i < nX; ++i) {buffer[i].rank = -1; buffer[i].index = -1;}
1357:     PetscSFReduceBegin(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1358:     PetscSFReduceEnd(sfXA, MPIU_2INT, owners, buffer, MPI_MAXLOC);
1359:     PetscSFDestroy(&sfXA);
1360:     PetscFree(owners);
1361:     for (i = 0, nleaves = 0; i < nX; ++i) if (buffer[i].rank >= 0) nleaves++;
1362:     PetscMalloc1(nleaves, &ilocal);
1363:     PetscMalloc1(nleaves, &iremote);
1364:     for (i = 0, nleaves = 0; i < nX; ++i) {
1365:       if (buffer[i].rank >= 0) {
1366:         ilocal[nleaves] = i;
1367:         iremote[nleaves].rank = buffer[i].rank;
1368:         iremote[nleaves].index = buffer[i].index;
1369:         nleaves++;
1370:       }
1371:     }
1372:     PetscSFCreate(comm, &sfAX);
1373:     PetscSFSetFromOptions(sfAX);
1374:     PetscSFSetGraph(sfAX, n, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
1375:     /* Fix PetscSFCompose() and replace the code-block below with:  */
1376:     /* PetscSFCompose(sfAX, sfXB, &sfAB);      */
1377:     /* which currently causes segmentation fault due to sparse map. */
1378:     {
1379:       PetscInt     npoints;
1380:       PetscInt     mleaves;
1381:       PetscInt    *jlocal;
1382:       PetscSFNode *jremote;

1384:       PetscSFGetGraph(sfXB, NULL, &npoints, NULL, NULL);
1385:       PetscMalloc1(npoints, &owners);
1386:       for (i = 0; i < npoints; ++i) {owners[i].rank = -1; owners[i].index = -1;}
1387:       PetscSFBcastBegin(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1388:       PetscSFBcastEnd(sfXB, MPIU_2INT, buffer, owners, MPI_REPLACE);
1389:       for (i = 0, mleaves = 0; i < npoints; ++i) if (owners[i].rank >= 0) mleaves++;
1390:       PetscMalloc1(mleaves, &jlocal);
1391:       PetscMalloc1(mleaves, &jremote);
1392:       for (i = 0, mleaves = 0; i < npoints; ++i) {
1393:         if (owners[i].rank >= 0) {
1394:           jlocal[mleaves] = i;
1395:           jremote[mleaves].rank = owners[i].rank;
1396:           jremote[mleaves].index = owners[i].index;
1397:           mleaves++;
1398:         }
1399:       }
1400:       PetscSFCreate(comm, &sfAB);
1401:       PetscSFSetFromOptions(sfAB);
1402:       PetscSFSetGraph(sfAB, n, mleaves, jlocal, PETSC_OWN_POINTER, jremote, PETSC_OWN_POINTER);
1403:       PetscFree(owners);
1404:     }
1405:     PetscFree(buffer);
1406:     PetscSFDestroy(&sfAX);
1407:   }
1408:   PetscViewerHDF5PopGroup(viewer);
1409:   PetscViewerHDF5PopGroup(viewer);
1410:   PetscViewerHDF5PopGroup(viewer);
1411:   PetscViewerHDF5PopGroup(viewer);
1412:   /* Create plex section (sectionB) */
1413:   DMGetLocalSection(sectiondm, &sectionB);
1414:   if (lsf || gsf) {
1415:     PetscLayout  layout;
1416:     PetscInt     M, m;
1417:     PetscInt    *offsetsA;
1418:     PetscBool    includesConstraintsA;

1420:     PetscSFDistributeSection(sfAB, sectionA, &offsetsA, sectionB);
1421:     PetscSectionGetIncludesConstraints(sectionA, &includesConstraintsA);
1422:     if (includesConstraintsA) {PetscSectionGetStorageSize(sectionA, &m);}
1423:     else {PetscSectionGetConstrainedStorageSize(sectionA, &m);}
1424:     MPI_Allreduce(&m, &M, 1, MPIU_INT, MPI_SUM, comm);
1425:     PetscLayoutCreate(comm, &layout);
1426:     PetscLayoutSetSize(layout, M);
1427:     PetscLayoutSetUp(layout);
1428:     if (lsf) {
1429:       PetscSF lsfABdata;

1431:       DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, sectionB, &lsfABdata);
1432:       *lsf = lsfABdata;
1433:     }
1434:     if (gsf) {
1435:       PetscSection  gsectionB, gsectionB1;
1436:       PetscBool     includesConstraintsB;
1437:       PetscSF       gsfABdata, pointsf;

1439:       DMGetGlobalSection(sectiondm, &gsectionB1);
1440:       PetscSectionGetIncludesConstraints(gsectionB1, &includesConstraintsB);
1441:       DMGetPointSF(sectiondm, &pointsf);
1442:       PetscSectionCreateGlobalSection(sectionB, pointsf, includesConstraintsB, PETSC_TRUE, &gsectionB);
1443:       DMPlexSectionLoad_HDF5_Internal_CreateDataSF(sectionA, layout, offsetsA, gsectionB, &gsfABdata);
1444:       PetscSectionDestroy(&gsectionB);
1445:       *gsf = gsfABdata;
1446:     }
1447:     PetscLayoutDestroy(&layout);
1448:     PetscFree(offsetsA);
1449:   } else {
1450:     PetscSFDistributeSection(sfAB, sectionA, NULL, sectionB);
1451:   }
1452:   PetscSFDestroy(&sfAB);
1453:   PetscSectionDestroy(&sectionA);
1454:   return(0);
1455: }

1457: PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM dm, PetscViewer viewer, DM sectiondm, PetscSF sf, Vec vec)
1458: {
1459:   MPI_Comm           comm;
1460:   const char        *topologydm_name;
1461:   const char        *sectiondm_name;
1462:   const char        *vec_name;
1463:   Vec                vecA;
1464:   PetscInt           mA, m, bs;
1465:   const PetscInt    *ilocal;
1466:   const PetscScalar *src;
1467:   PetscScalar       *dest;
1468:   PetscErrorCode     ierr;

1471:   PetscObjectGetComm((PetscObject)dm, &comm);
1472:   PetscObjectGetName((PetscObject)dm, &topologydm_name);
1473:   PetscObjectGetName((PetscObject)sectiondm, &sectiondm_name);
1474:   PetscObjectGetName((PetscObject)vec, &vec_name);
1475:   PetscViewerHDF5PushGroup(viewer, "topologies");
1476:   PetscViewerHDF5PushGroup(viewer, topologydm_name);
1477:   PetscViewerHDF5PushGroup(viewer, "dms");
1478:   PetscViewerHDF5PushGroup(viewer, sectiondm_name);
1479:   PetscViewerHDF5PushGroup(viewer, "vecs");
1480:   PetscViewerHDF5PushGroup(viewer, vec_name);
1481:   VecCreate(comm, &vecA);
1482:   PetscObjectSetName((PetscObject)vecA, vec_name);
1483:   PetscSFGetGraph(sf, &mA, &m, &ilocal, NULL);
1484:   /* Check consistency */
1485:   {
1486:     PetscSF   pointsf, pointsf1;
1487:     PetscInt  m1, i, j;

1489:     DMGetPointSF(dm, &pointsf);
1490:     DMGetPointSF(sectiondm, &pointsf1);
1491:     if (pointsf1 != pointsf) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatching point SFs for dm and sectiondm");
1492: #if defined(PETSC_USE_DEBUG)
1493:     {
1494:       PetscInt  MA, MA1;

1496:       MPIU_Allreduce(&mA, &MA, 1, MPIU_INT, MPI_SUM, comm);
1497:       PetscViewerHDF5ReadSizes(viewer, vec_name, NULL, &MA1);
1498:       if (MA1 != MA) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total SF root size (%D) != On-disk vector data size (%D)", MA, MA1);
1499:     }
1500: #endif
1501:     VecGetLocalSize(vec, &m1);
1502:     if (m1 < m) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Target vector size (%D) < SF leaf size (%D)", m1, m);
1503:     for (i = 0; i < m; ++i) {
1504:       j = ilocal ? ilocal[i] : i;
1505:       if (j < 0 || j >= m1) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Leaf's %D-th index, %D, not in [%D, %D)", i, j, 0, m1);
1506:     }
1507:   }
1508:   VecSetSizes(vecA, mA, PETSC_DECIDE);
1509:   VecLoad(vecA, viewer);
1510:   VecGetArrayRead(vecA, &src);
1511:   VecGetArray(vec, &dest);
1512:   PetscSFBcastBegin(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1513:   PetscSFBcastEnd(sf, MPIU_SCALAR, src, dest, MPI_REPLACE);
1514:   VecRestoreArray(vec, &dest);
1515:   VecRestoreArrayRead(vecA, &src);
1516:   VecDestroy(&vecA);
1517:   PetscViewerHDF5ReadAttribute(viewer, NULL, "blockSize", PETSC_INT, NULL, (void *) &bs);
1518:   VecSetBlockSize(vec, bs);
1519:   PetscViewerHDF5PopGroup(viewer);
1520:   PetscViewerHDF5PopGroup(viewer);
1521:   PetscViewerHDF5PopGroup(viewer);
1522:   PetscViewerHDF5PopGroup(viewer);
1523:   PetscViewerHDF5PopGroup(viewer);
1524:   PetscViewerHDF5PopGroup(viewer);
1525:   return(0);
1526: }
1527: #endif