Actual source code: plexhdf5.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petsc/private/isimpl.h>
  3: #include <petsc/private/vecimpl.h>
  4: #include <petscviewerhdf5.h>

  6: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

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

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

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

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

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

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

 76: PetscErrorCode VecView_Plex_Local_HDF5(Vec v, PetscViewer viewer)
 77: {
 78:   DM                      dm;
 79:   DM                      dmBC;
 80:   PetscSection            section, sectionGlobal;
 81:   Vec                     gv;
 82:   const char             *name;
 83:   PetscViewerVTKFieldType ft;
 84:   PetscViewerFormat       format;
 85:   PetscInt                seqnum;
 86:   PetscReal               seqval;
 87:   PetscBool               isseq;
 88:   PetscErrorCode          ierr;

 91:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
 92:   VecGetDM(v, &dm);
 93:   DMGetDefaultSection(dm, &section);
 94:   DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
 95:   PetscViewerHDF5SetTimestep(viewer, seqnum);
 96:   DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
 97:   PetscViewerGetFormat(viewer, &format);
 98:   DMGetOutputDM(dm, &dmBC);
 99:   DMGetDefaultGlobalSection(dmBC, &sectionGlobal);
100:   DMGetGlobalVector(dmBC, &gv);
101:   PetscObjectGetName((PetscObject) v, &name);
102:   PetscObjectSetName((PetscObject) gv, name);
103:   DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
104:   DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
105:   PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
106:   if (isseq) {VecView_Seq(gv, viewer);}
107:   else       {VecView_MPI(gv, viewer);}
108:   if (format == PETSC_VIEWER_HDF5_VIZ) {
109:     /* Output visualization representation */
110:     PetscInt numFields, f;

112:     PetscSectionGetNumFields(section, &numFields);
113:     for (f = 0; f < numFields; ++f) {
114:       Vec         subv;
115:       IS          is;
116:       const char *fname, *fgroup;
117:       char        subname[PETSC_MAX_PATH_LEN];
118:       char        group[PETSC_MAX_PATH_LEN];
119:       PetscInt    pStart, pEnd;
120:       PetscBool   flag;

122:       DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
123:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
124:       PetscSectionGetFieldName(section, f, &fname);
125:       if (!fname) continue;
126:       PetscViewerHDF5PushGroup(viewer, fgroup);
127:       PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
128:       PetscStrcpy(subname, name);
129:       PetscStrcat(subname, "_");
130:       PetscStrcat(subname, fname);
131:       PetscObjectSetName((PetscObject) subv, subname);
132:       if (isseq) {VecView_Seq(subv, viewer);}
133:       else       {VecView_MPI(subv, viewer);}
134:       PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
135:       PetscViewerHDF5PopGroup(viewer);
136:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s/%s", fgroup, subname);
137:       PetscViewerHDF5HasAttribute(viewer, group, "vector_field_type", &flag);
138:       if (!flag) {
139:         if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
140:           PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "vector");
141:         } else {
142:           PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "scalar");
143:         }
144:       }
145:     }
146:   }
147:   DMRestoreGlobalVector(dmBC, &gv);
148:   return(0);
149: }

153: PetscErrorCode VecView_Plex_HDF5(Vec v, PetscViewer viewer)
154: {
155:   DM             dm;
156:   Vec            locv;
157:   const char    *name;

161:   VecGetDM(v, &dm);
162:   DMGetLocalVector(dm, &locv);
163:   PetscObjectGetName((PetscObject) v, &name);
164:   PetscObjectSetName((PetscObject) locv, name);
165:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
166:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
167:   PetscViewerHDF5PushGroup(viewer, "/fields");
168:   PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
169:   VecView_Plex_Local(locv, viewer);
170:   PetscViewerPopFormat(viewer);
171:   PetscViewerHDF5PopGroup(viewer);
172:   DMRestoreLocalVector(dm, &locv);
173:   return(0);
174: }

178: PetscErrorCode VecView_Plex_HDF5_Native(Vec v, PetscViewer viewer)
179: {
180:   PetscBool      isseq;

184:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
185:   PetscViewerHDF5PushGroup(viewer, "/fields");
186:   if (isseq) {VecView_Seq(v, viewer);}
187:   else       {VecView_MPI(v, viewer);}
188:   PetscViewerHDF5PopGroup(viewer);
189:   return(0);
190: }

194: PetscErrorCode VecLoad_Plex_HDF5(Vec v, PetscViewer viewer)
195: {
196:   DM             dm;
197:   Vec            locv;
198:   const char    *name;
199:   PetscInt       seqnum;

203:   VecGetDM(v, &dm);
204:   DMGetLocalVector(dm, &locv);
205:   PetscObjectGetName((PetscObject) v, &name);
206:   PetscObjectSetName((PetscObject) locv, name);
207:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
208:   PetscViewerHDF5SetTimestep(viewer, seqnum);
209:   PetscViewerHDF5PushGroup(viewer, "/fields");
210:   VecLoad_Plex_Local(locv, viewer);
211:   PetscViewerHDF5PopGroup(viewer);
212:   DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
213:   DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
214:   DMRestoreLocalVector(dm, &locv);
215:   return(0);
216: }

220: PetscErrorCode VecLoad_Plex_HDF5_Native(Vec v, PetscViewer viewer)
221: {
222:   DM             dm;
223:   PetscInt       seqnum;

227:   VecGetDM(v, &dm);
228:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
229:   PetscViewerHDF5SetTimestep(viewer, seqnum);
230:   PetscViewerHDF5PushGroup(viewer, "/fields");
231:   VecLoad_Default(v, viewer);
232:   PetscViewerHDF5PopGroup(viewer);
233:   return(0);
234: }

238: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
239: {
240:   IS              orderIS, conesIS, cellsIS, orntsIS;
241:   const PetscInt *gpoint;
242:   PetscInt       *order, *sizes, *cones, *ornts;
243:   PetscInt        dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
244:   PetscErrorCode  ierr;

247:   ISGetIndices(globalPointNumbers, &gpoint);
248:   DMGetDimension(dm, &dim);
249:   DMPlexGetChart(dm, &pStart, &pEnd);
250:   for (p = pStart; p < pEnd; ++p) {
251:     if (gpoint[p] >= 0) {
252:       PetscInt coneSize;

254:       DMPlexGetConeSize(dm, p, &coneSize);
255:       conesSize += 1;
256:       cellsSize += coneSize;
257:     }
258:   }
259:   PetscMalloc1(conesSize, &order);
260:   PetscMalloc1(conesSize, &sizes);
261:   PetscMalloc1(cellsSize, &cones);
262:   PetscMalloc1(cellsSize, &ornts);
263:   for (p = pStart; p < pEnd; ++p) {
264:     if (gpoint[p] >= 0) {
265:       const PetscInt *cone, *ornt;
266:       PetscInt        coneSize, cp;

268:       DMPlexGetConeSize(dm, p, &coneSize);
269:       DMPlexGetCone(dm, p, &cone);
270:       DMPlexGetConeOrientation(dm, p, &ornt);
271:       order[s]   = gpoint[p];
272:       sizes[s++] = coneSize;
273:       for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
274:     }
275:   }
276:   if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %d != %d", s, conesSize);
277:   if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %d != %d", c, cellsSize);
278:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
279:   PetscObjectSetName((PetscObject) orderIS, "order");
280:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
281:   PetscObjectSetName((PetscObject) conesIS, "cones");
282:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
283:   PetscObjectSetName((PetscObject) cellsIS, "cells");
284:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
285:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
286:   PetscViewerHDF5PushGroup(viewer, "/topology");
287:   ISView(orderIS, viewer);
288:   ISView(conesIS, viewer);
289:   ISView(cellsIS, viewer);
290:   ISView(orntsIS, viewer);
291:   PetscViewerHDF5PopGroup(viewer);
292:   ISDestroy(&orderIS);
293:   ISDestroy(&conesIS);
294:   ISDestroy(&cellsIS);
295:   ISDestroy(&orntsIS);
296:   ISRestoreIndices(globalPointNumbers, &gpoint);

298:   PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
299:   return(0);
300: }

304: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, DMLabel label, PetscInt labelId, PetscViewer viewer)
305: {
306:   IS              cellIS, globalVertexNumbers;
307:   const PetscInt *gvertex;
308:   PetscInt       *vertices;
309:   PetscInt        dim, depth, vStart, vEnd, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
310:   hid_t           fileId, groupId;
311:   PetscErrorCode  ierr;

314:   DMGetDimension(dm, &dim);
315:   DMPlexGetDepth(dm, &depth);
316:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
317:   DMPlexGetVTKCellHeight(dm, &cellHeight);
318:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
319:   DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, PETSC_NULL);
320:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
321:   if (depth == 1) return(0);
322:   for (cell = cStart; cell < cEnd; ++cell) {
323:     PetscInt *closure = NULL;
324:     PetscInt  closureSize, v, Nc = 0;

326:     if (label) {
327:       PetscInt value;
328:       DMLabelGetValue(label, cell, &value);
329:       if (value == labelId) continue;
330:     }
331:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
332:     for (v = 0; v < closureSize*2; v += 2) {
333:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
334:     }
335:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
336:     conesSize += Nc;
337:     if (!numCornersLocal)           numCornersLocal = Nc;
338:     else if (numCornersLocal != Nc) numCornersLocal = 1;
339:   }
340:   MPIU_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
341:   if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");

343:   PetscViewerHDF5PushGroup(viewer, "/viz");
344:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
345:   PetscStackCallHDF5(H5Gclose,(groupId));
346:   PetscViewerHDF5PopGroup(viewer);

348:   DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
349:   ISGetIndices(globalVertexNumbers, &gvertex);
350:   PetscMalloc1(conesSize, &vertices);
351:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
352:     PetscInt *closure = NULL;
353:     PetscInt  closureSize, Nc = 0, p;

355:     if (label) {
356:       PetscInt value;
357:       DMLabelGetValue(label, cell, &value);
358:       if (value == labelId) continue;
359:     }
360:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
361:     for (p = 0; p < closureSize*2; p += 2) {
362:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
363:         closure[Nc++] = closure[p];
364:       }
365:     }
366:     DMPlexInvertCell_Internal(dim, Nc, closure);
367:     for (p = 0; p < Nc; ++p) {
368:       const PetscInt gv = gvertex[closure[p] - vStart];
369:       vertices[v++] = gv < 0 ? -(gv+1) : gv;
370:     }
371:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
372:   }
373:   if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
374:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);
375:   PetscLayoutSetBlockSize(cellIS->map, numCorners);
376:   PetscObjectSetName((PetscObject) cellIS, "cells");
377:   PetscViewerHDF5PushGroup(viewer, "/viz/topology");
378:   ISView(cellIS, viewer);
379: #if 0
380:   if (numCorners == 1) {
381:     ISView(coneIS, viewer);
382:   } else {
383:     PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
384:   }
385:   ISDestroy(&coneIS);
386: #else
387:   PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
388: #endif
389:   ISDestroy(&cellIS);

391:   PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
392:   return(0);
393: }

397: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
398: {
399:   DM             cdm;
400:   Vec            coordinates, newcoords;
401:   PetscReal      lengthScale;
402:   PetscInt       m, M, bs;

406:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
407:   DMGetCoordinateDM(dm, &cdm);
408:   DMGetCoordinates(dm, &coordinates);
409:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
410:   PetscObjectSetName((PetscObject) newcoords, "vertices");
411:   VecGetSize(coordinates, &M);
412:   VecGetLocalSize(coordinates, &m);
413:   VecSetSizes(newcoords, m, M);
414:   VecGetBlockSize(coordinates, &bs);
415:   VecSetBlockSize(newcoords, bs);
416:   VecSetType(newcoords,VECSTANDARD);
417:   VecCopy(coordinates, newcoords);
418:   VecScale(newcoords, lengthScale);
419:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
420:   PetscViewerHDF5PushGroup(viewer, "/geometry");
421:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
422:   VecView(newcoords, viewer);
423:   PetscViewerPopFormat(viewer);
424:   PetscViewerHDF5PopGroup(viewer);
425:   VecDestroy(&newcoords);
426:   return(0);
427: }

431: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
432: {
433:   Vec              coordinates, newcoords;
434:   PetscSection     cSection;
435:   PetscScalar     *coords, *ncoords;
436:   const PetscReal *L;
437:   const DMBoundaryType *bd;
438:   PetscReal        lengthScale;
439:   PetscInt         vStart, vEnd, v, bs, coordSize, dof, off, d;
440:   hid_t            fileId, groupId;
441:   PetscErrorCode   ierr;

444:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
445:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
446:   DMGetCoordinatesLocal(dm, &coordinates);
447:   VecGetBlockSize(coordinates, &bs);
448:   VecGetLocalSize(coordinates, &coordSize);
449:   if (coordSize == (vEnd - vStart)*bs) return(0);
450:   DMGetPeriodicity(dm, NULL, &L, &bd);
451:   DMGetCoordinateSection(dm, &cSection);
452:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
453:   coordSize = 0;
454:   for (v = vStart; v < vEnd; ++v) {
455:     PetscSectionGetDof(cSection, v, &dof);
456:     if (L && dof == 2) coordSize += dof+1;
457:     else               coordSize += dof;
458:   }
459:   if (L && bs == 2) {VecSetBlockSize(newcoords, bs+1);}
460:   else              {VecSetBlockSize(newcoords, bs);}
461:   VecSetSizes(newcoords, coordSize, PETSC_DETERMINE);
462:   VecSetType(newcoords,VECSTANDARD);
463:   VecGetArray(coordinates, &coords);
464:   VecGetArray(newcoords,   &ncoords);
465:   coordSize = 0;
466:   for (v = vStart; v < vEnd; ++v) {
467:     PetscSectionGetDof(cSection, v, &dof);
468:     PetscSectionGetOffset(cSection, v, &off);
469:     if (L && dof == 2) {
470:       /* Need to do torus */
471:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot embed doubly periodic domain yet");}
472:       else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
473:         /* X-periodic */
474:         ncoords[coordSize++] = -cos(2.0*PETSC_PI*coords[off+0]/L[0])*(L[0]/(2.0*PETSC_PI));
475:         ncoords[coordSize++] = coords[off+1];
476:         ncoords[coordSize++] = sin(2.0*PETSC_PI*coords[off+0]/L[0])*(L[0]/(2.0*PETSC_PI));
477:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
478:         /* Y-periodic */
479:         ncoords[coordSize++] = coords[off+0];
480:         ncoords[coordSize++] = sin(2.0*PETSC_PI*coords[off+1]/L[1])*(L[1]/(2.0*PETSC_PI));
481:         ncoords[coordSize++] = -cos(2.0*PETSC_PI*coords[off+1]/L[1])*(L[1]/(2.0*PETSC_PI));
482:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
483:     } else {
484:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
485:     }
486:   }
487:   VecRestoreArray(coordinates, &coords);
488:   VecRestoreArray(newcoords,   &ncoords);
489:   PetscObjectSetName((PetscObject) newcoords, "vertices");
490:   VecScale(newcoords, lengthScale);
491:   PetscViewerHDF5PushGroup(viewer, "/viz");
492:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
493:   PetscStackCallHDF5(H5Gclose,(groupId));
494:   PetscViewerHDF5PopGroup(viewer);
495:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
496:   VecView(newcoords, viewer);
497:   PetscViewerHDF5PopGroup(viewer);
498:   VecDestroy(&newcoords);
499:   return(0);
500: }

504: /* We only write cells and vertices. Does this screw up parallel reading? */
505: PetscErrorCode DMPlexView_HDF5(DM dm, PetscViewer viewer)
506: {
507:   DMLabel           label   = NULL;
508:   PetscInt          labelId = 0;
509:   IS                globalPointNumbers;
510:   const PetscInt   *gpoint;
511:   PetscInt          numLabels, l;
512:   hid_t             fileId, groupId;
513:   PetscViewerFormat format;
514:   PetscErrorCode    ierr;

517:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
518:   PetscViewerGetFormat(viewer, &format);
519:   DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
520:   if (format == PETSC_VIEWER_HDF5_VIZ) {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
521:   DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);
522:   if (format == PETSC_VIEWER_HDF5_VIZ) {DMPlexWriteTopology_Vertices_HDF5_Static(dm, label, labelId, viewer);}
523:   /* Write Labels*/
524:   ISGetIndices(globalPointNumbers, &gpoint);
525:   PetscViewerHDF5PushGroup(viewer, "/labels");
526:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
527:   if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
528:   PetscViewerHDF5PopGroup(viewer);
529:   DMGetNumLabels(dm, &numLabels);
530:   for (l = 0; l < numLabels; ++l) {
531:     DMLabel         label;
532:     const char     *name;
533:     IS              valueIS, pvalueIS, globalValueIS;
534:     const PetscInt *values;
535:     PetscInt        numValues, v;
536:     PetscBool       isDepth, output;
537:     char            group[PETSC_MAX_PATH_LEN];

539:     DMGetLabelName(dm, l, &name);
540:     DMGetLabelOutput(dm, name, &output);
541:     PetscStrncmp(name, "depth", 10, &isDepth);
542:     if (isDepth || !output) continue;
543:     DMGetLabel(dm, name, &label);
544:     PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
545:     PetscViewerHDF5PushGroup(viewer, group);
546:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
547:     if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
548:     PetscViewerHDF5PopGroup(viewer);
549:     DMLabelGetValueIS(label, &valueIS);
550:     /* Must copy to a new IS on the global comm */
551:     ISGetLocalSize(valueIS, &numValues);
552:     ISGetIndices(valueIS, &values);
553:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
554:     ISRestoreIndices(valueIS, &values);
555:     ISAllGather(pvalueIS, &globalValueIS);
556:     ISDestroy(&pvalueIS);
557:     ISSortRemoveDups(globalValueIS);
558:     ISGetLocalSize(globalValueIS, &numValues);
559:     ISGetIndices(globalValueIS, &values);
560:     for (v = 0; v < numValues; ++v) {
561:       IS              stratumIS, globalStratumIS;
562:       const PetscInt *spoints = NULL;
563:       PetscInt       *gspoints, n = 0, gn, p;
564:       const char     *iname = "indices";

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

569:       if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
570:       if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
571:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
572:       PetscMalloc1(gn,&gspoints);
573:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
574:       if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
575:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
576:       if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
577:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

579:       PetscViewerHDF5PushGroup(viewer, group);
580:       ISView(globalStratumIS, viewer);
581:       PetscViewerHDF5PopGroup(viewer);
582:       ISDestroy(&globalStratumIS);
583:       ISDestroy(&stratumIS);
584:     }
585:     ISRestoreIndices(globalValueIS, &values);
586:     ISDestroy(&globalValueIS);
587:     ISDestroy(&valueIS);
588:   }
589:   ISRestoreIndices(globalPointNumbers, &gpoint);
590:   ISDestroy(&globalPointNumbers);
591:   return(0);
592: }

594: typedef struct {
595:   PetscMPIInt rank;
596:   DM          dm;
597:   PetscViewer viewer;
598:   DMLabel     label;
599: } LabelCtx;

601: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
602: {
603:   PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
604:   DMLabel         label  = ((LabelCtx *) op_data)->label;
605:   IS              stratumIS;
606:   const PetscInt *ind;
607:   PetscInt        value, N, i;
608:   const char     *lname;
609:   char            group[PETSC_MAX_PATH_LEN];
610:   PetscErrorCode  ierr;

612:   PetscOptionsStringToInt(name, &value);
613:   ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
614:   PetscObjectSetName((PetscObject) stratumIS, "indices");
615:   DMLabelGetName(label, &lname);
616:   PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
617:   PetscViewerHDF5PushGroup(viewer, group);
618:   {
619:     /* Force serial load */
620:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
621:     PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
622:     PetscLayoutSetSize(stratumIS->map, N);
623:   }
624:   ISLoad(stratumIS, viewer);
625:   PetscViewerHDF5PopGroup(viewer);
626:   ISGetLocalSize(stratumIS, &N);
627:   ISGetIndices(stratumIS, &ind);
628:   for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
629:   ISRestoreIndices(stratumIS, &ind);
630:   ISDestroy(&stratumIS);
631:   return 0;
632: }

634: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
635: {
636:   DM             dm  = ((LabelCtx *) op_data)->dm;
637:   hsize_t        idx = 0;
639:   herr_t         err;

641:   DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
642:   DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
643:   PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
644:   return err;
645: }

649: /* The first version will read everything onto proc 0, letting the user distribute
650:    The next will create a naive partition, and then rebalance after reading
651: */
652: PetscErrorCode DMPlexLoad_HDF5(DM dm, PetscViewer viewer)
653: {
654:   LabelCtx        ctx;
655:   PetscSection    coordSection;
656:   Vec             coordinates;
657:   IS              orderIS, conesIS, cellsIS, orntsIS;
658:   const PetscInt *order, *cones, *cells, *ornts;
659:   PetscReal       lengthScale;
660:   PetscInt       *cone, *ornt;
661:   PetscInt        dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
662:   hid_t           fileId, groupId;
663:   hsize_t         idx = 0;
664:   PetscMPIInt     rank;
665:   PetscErrorCode  ierr;

668:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
669:   /* Read toplogy */
670:   PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
671:   DMSetDimension(dm, dim);
672:   PetscViewerHDF5PushGroup(viewer, "/topology");

674:   ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
675:   PetscObjectSetName((PetscObject) orderIS, "order");
676:   ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
677:   PetscObjectSetName((PetscObject) conesIS, "cones");
678:   ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
679:   PetscObjectSetName((PetscObject) cellsIS, "cells");
680:   ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
681:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
682:   {
683:     /* Force serial load */
684:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
685:     PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
686:     PetscLayoutSetSize(orderIS->map, pEnd);
687:     PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
688:     PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
689:     PetscLayoutSetSize(conesIS->map, pEnd);
690:     pEnd = !rank ? pEnd : 0;
691:     PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
692:     PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
693:     PetscLayoutSetSize(cellsIS->map, N);
694:     PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
695:     PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
696:     PetscLayoutSetSize(orntsIS->map, N);
697:   }
698:   ISLoad(orderIS, viewer);
699:   ISLoad(conesIS, viewer);
700:   ISLoad(cellsIS, viewer);
701:   ISLoad(orntsIS, viewer);
702:   PetscViewerHDF5PopGroup(viewer);
703:   /* Read geometry */
704:   PetscViewerHDF5PushGroup(viewer, "/geometry");
705:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
706:   PetscObjectSetName((PetscObject) coordinates, "vertices");
707:   {
708:     /* Force serial load */
709:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
710:     VecSetSizes(coordinates, !rank ? N : 0, N);
711:     VecSetBlockSize(coordinates, spatialDim);
712:   }
713:   VecLoad(coordinates, viewer);
714:   PetscViewerHDF5PopGroup(viewer);
715:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
716:   VecScale(coordinates, 1.0/lengthScale);
717:   VecGetLocalSize(coordinates, &numVertices);
718:   VecGetBlockSize(coordinates, &spatialDim);
719:   numVertices /= spatialDim;
720:   /* Create Plex */
721:   DMPlexSetChart(dm, 0, pEnd);
722:   ISGetIndices(orderIS, &order);
723:   ISGetIndices(conesIS, &cones);
724:   for (p = 0; p < pEnd; ++p) {
725:     DMPlexSetConeSize(dm, order[p], cones[p]);
726:     maxConeSize = PetscMax(maxConeSize, cones[p]);
727:   }
728:   DMSetUp(dm);
729:   ISGetIndices(cellsIS, &cells);
730:   ISGetIndices(orntsIS, &ornts);
731:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
732:   for (p = 0, q = 0; p < pEnd; ++p) {
733:     for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
734:     DMPlexSetCone(dm, order[p], cone);
735:     DMPlexSetConeOrientation(dm, order[p], ornt);
736:   }
737:   PetscFree2(cone,ornt);
738:   ISRestoreIndices(orderIS, &order);
739:   ISRestoreIndices(conesIS, &cones);
740:   ISRestoreIndices(cellsIS, &cells);
741:   ISRestoreIndices(orntsIS, &ornts);
742:   ISDestroy(&orderIS);
743:   ISDestroy(&conesIS);
744:   ISDestroy(&cellsIS);
745:   ISDestroy(&orntsIS);
746:   DMPlexSymmetrize(dm);
747:   DMPlexStratify(dm);
748:   /* Create coordinates */
749:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
750:   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);
751:   DMGetCoordinateSection(dm, &coordSection);
752:   PetscSectionSetNumFields(coordSection, 1);
753:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
754:   PetscSectionSetChart(coordSection, vStart, vEnd);
755:   for (v = vStart; v < vEnd; ++v) {
756:     PetscSectionSetDof(coordSection, v, spatialDim);
757:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
758:   }
759:   PetscSectionSetUp(coordSection);
760:   DMSetCoordinates(dm, coordinates);
761:   VecDestroy(&coordinates);
762:   /* Read Labels*/
763:   ctx.rank   = rank;
764:   ctx.dm     = dm;
765:   ctx.viewer = viewer;
766:   PetscViewerHDF5PushGroup(viewer, "/labels");
767:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
768:   PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
769:   PetscStackCallHDF5(H5Gclose,(groupId));
770:   PetscViewerHDF5PopGroup(viewer);
771:   return(0);
772: }
773: #endif