Actual source code: plexhdf5.c

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

  6: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);

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

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

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

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

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

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

 70: PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec v, PetscViewer viewer)
 71: {
 72:   DM                      dm;
 73:   DM                      dmBC;
 74:   PetscSection            section, sectionGlobal;
 75:   Vec                     gv;
 76:   const char             *name;
 77:   PetscViewerVTKFieldType ft;
 78:   PetscViewerFormat       format;
 79:   PetscInt                seqnum;
 80:   PetscReal               seqval;
 81:   PetscBool               isseq;
 82:   PetscErrorCode          ierr;

 85:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
 86:   VecGetDM(v, &dm);
 87:   DMGetSection(dm, &section);
 88:   DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
 89:   PetscViewerHDF5SetTimestep(viewer, seqnum);
 90:   DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
 91:   PetscViewerGetFormat(viewer, &format);
 92:   DMGetOutputDM(dm, &dmBC);
 93:   DMGetGlobalSection(dmBC, &sectionGlobal);
 94:   DMGetGlobalVector(dmBC, &gv);
 95:   PetscObjectGetName((PetscObject) v, &name);
 96:   PetscObjectSetName((PetscObject) gv, name);
 97:   DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
 98:   DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
 99:   PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
100:   if (isseq) {VecView_Seq(gv, viewer);}
101:   else       {VecView_MPI(gv, viewer);}
102:   if (format == PETSC_VIEWER_HDF5_VIZ) {
103:     /* Output visualization representation */
104:     PetscInt numFields, f;
105:     DMLabel  cutLabel;

107:     PetscSectionGetNumFields(section, &numFields);
108:     DMGetLabel(dm, "periodic_cut", &cutLabel);
109:     for (f = 0; f < numFields; ++f) {
110:       Vec         subv;
111:       IS          is;
112:       const char *fname, *fgroup;
113:       char        subname[PETSC_MAX_PATH_LEN];
114:       char        group[PETSC_MAX_PATH_LEN];
115:       PetscInt    pStart, pEnd;
116:       PetscBool   flag;

118:       DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
119:       fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
120:       PetscSectionGetFieldName(section, f, &fname);
121:       if (!fname) continue;
122:       PetscViewerHDF5PushGroup(viewer, fgroup);
123:       if (cutLabel) {
124:         const PetscScalar *ga;
125:         PetscScalar       *suba;
126:         PetscInt           Nc, gstart, subSize = 0, extSize = 0, subOff = 0, newOff = 0, p;

128:         PetscSectionGetFieldComponents(section, f, &Nc);
129:         for (p = pStart; p < pEnd; ++p) {
130:           PetscInt gdof, fdof = 0, val;

132:           PetscSectionGetDof(sectionGlobal, p, &gdof);
133:           if (gdof > 0) {PetscSectionGetFieldDof(section, p, f, &fdof);}
134:           subSize += fdof;
135:           DMLabelGetValue(cutLabel, p, &val);
136:           if (val == 1) extSize += fdof;
137:         }
138:         VecCreate(PetscObjectComm((PetscObject) gv), &subv);
139:         VecSetSizes(subv, subSize+extSize, PETSC_DETERMINE);
140:         VecSetBlockSize(subv, Nc);
141:         VecSetType(subv, VECSTANDARD);
142:         VecGetOwnershipRange(gv, &gstart, NULL);
143:         VecGetArrayRead(gv, &ga);
144:         VecGetArray(subv, &suba);
145:         for (p = pStart; p < pEnd; ++p) {
146:           PetscInt gdof, goff, val;

148:           PetscSectionGetDof(sectionGlobal, p, &gdof);
149:           if (gdof > 0) {
150:             PetscInt fdof, fc, f2, poff = 0;

152:             PetscSectionGetOffset(sectionGlobal, p, &goff);
153:             /* Can get rid of this loop by storing field information in the global section */
154:             for (f2 = 0; f2 < f; ++f2) {
155:               PetscSectionGetFieldDof(section, p, f2, &fdof);
156:               poff += fdof;
157:             }
158:             PetscSectionGetFieldDof(section, p, f, &fdof);
159:             for (fc = 0; fc < fdof; ++fc, ++subOff) suba[subOff] = ga[goff+poff+fc - gstart];
160:             DMLabelGetValue(cutLabel, p, &val);
161:             if (val == 1) {
162:               for (fc = 0; fc < fdof; ++fc, ++newOff) suba[subSize+newOff] = ga[goff+poff+fc - gstart];
163:             }
164:           }
165:         }
166:         VecRestoreArrayRead(gv, &ga);
167:         VecRestoreArray(subv, &suba);
168:       } else {
169:         PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
170:       }
171:       PetscStrncpy(subname, name,sizeof(subname));
172:       PetscStrlcat(subname, "_",sizeof(subname));
173:       PetscStrlcat(subname, fname,sizeof(subname));
174:       PetscObjectSetName((PetscObject) subv, subname);
175:       if (isseq) {VecView_Seq(subv, viewer);}
176:       else       {VecView_MPI(subv, viewer);}
177:       if (cutLabel) {
178:         VecDestroy(&subv);
179:       } else {
180:         PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
181:       }
182:       PetscViewerHDF5PopGroup(viewer);
183:       PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s/%s", fgroup, subname);
184:       PetscViewerHDF5HasAttribute(viewer, group, "vector_field_type", &flag);
185:       if (!flag) {
186:         if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
187:           PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "vector");
188:         } else {
189:           PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "scalar");
190:         }
191:       }
192:     }
193:   }
194:   DMRestoreGlobalVector(dmBC, &gv);
195:   return(0);
196: }

198: PetscErrorCode VecView_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
199: {
200:   DM             dm;
201:   Vec            locv;
202:   const char    *name;
203:   PetscReal      time;

207:   VecGetDM(v, &dm);
208:   DMGetLocalVector(dm, &locv);
209:   PetscObjectGetName((PetscObject) v, &name);
210:   PetscObjectSetName((PetscObject) locv, name);
211:   DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
212:   DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
213:   DMGetOutputSequenceNumber(dm, NULL, &time);
214:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locv, time, NULL, NULL, NULL);
215:   PetscViewerHDF5PushGroup(viewer, "/fields");
216:   PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
217:   VecView_Plex_Local_HDF5_Internal(locv, viewer);
218:   PetscViewerPopFormat(viewer);
219:   PetscViewerHDF5PopGroup(viewer);
220:   DMRestoreLocalVector(dm, &locv);
221:   return(0);
222: }

224: PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
225: {
226:   PetscBool      isseq;

230:   PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
231:   PetscViewerHDF5PushGroup(viewer, "/fields");
232:   if (isseq) {VecView_Seq(v, viewer);}
233:   else       {VecView_MPI(v, viewer);}
234:   PetscViewerHDF5PopGroup(viewer);
235:   return(0);
236: }

238: PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec v, PetscViewer viewer)
239: {
240:   DM             dm;
241:   Vec            locv;
242:   const char    *name;
243:   PetscInt       seqnum;

247:   VecGetDM(v, &dm);
248:   DMGetLocalVector(dm, &locv);
249:   PetscObjectGetName((PetscObject) v, &name);
250:   PetscObjectSetName((PetscObject) locv, name);
251:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
252:   PetscViewerHDF5SetTimestep(viewer, seqnum);
253:   PetscViewerHDF5PushGroup(viewer, "/fields");
254:   VecLoad_Plex_Local(locv, viewer);
255:   PetscViewerHDF5PopGroup(viewer);
256:   DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
257:   DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
258:   DMRestoreLocalVector(dm, &locv);
259:   return(0);
260: }

262: PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec v, PetscViewer viewer)
263: {
264:   DM             dm;
265:   PetscInt       seqnum;

269:   VecGetDM(v, &dm);
270:   DMGetOutputSequenceNumber(dm, &seqnum, NULL);
271:   PetscViewerHDF5SetTimestep(viewer, seqnum);
272:   PetscViewerHDF5PushGroup(viewer, "/fields");
273:   VecLoad_Default(v, viewer);
274:   PetscViewerHDF5PopGroup(viewer);
275:   return(0);
276: }

278: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
279: {
280:   IS              orderIS, conesIS, cellsIS, orntsIS;
281:   const PetscInt *gpoint;
282:   PetscInt       *order, *sizes, *cones, *ornts;
283:   PetscInt        dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
284:   PetscErrorCode  ierr;

287:   ISGetIndices(globalPointNumbers, &gpoint);
288:   DMGetDimension(dm, &dim);
289:   DMPlexGetChart(dm, &pStart, &pEnd);
290:   for (p = pStart; p < pEnd; ++p) {
291:     if (gpoint[p] >= 0) {
292:       PetscInt coneSize;

294:       DMPlexGetConeSize(dm, p, &coneSize);
295:       conesSize += 1;
296:       cellsSize += coneSize;
297:     }
298:   }
299:   PetscMalloc1(conesSize, &order);
300:   PetscMalloc1(conesSize, &sizes);
301:   PetscMalloc1(cellsSize, &cones);
302:   PetscMalloc1(cellsSize, &ornts);
303:   for (p = pStart; p < pEnd; ++p) {
304:     if (gpoint[p] >= 0) {
305:       const PetscInt *cone, *ornt;
306:       PetscInt        coneSize, cp;

308:       DMPlexGetConeSize(dm, p, &coneSize);
309:       DMPlexGetCone(dm, p, &cone);
310:       DMPlexGetConeOrientation(dm, p, &ornt);
311:       order[s]   = gpoint[p];
312:       sizes[s++] = coneSize;
313:       for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
314:     }
315:   }
316:   if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %d != %d", s, conesSize);
317:   if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %d != %d", c, cellsSize);
318:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
319:   PetscObjectSetName((PetscObject) orderIS, "order");
320:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
321:   PetscObjectSetName((PetscObject) conesIS, "cones");
322:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
323:   PetscObjectSetName((PetscObject) cellsIS, "cells");
324:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
325:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
326:   PetscViewerHDF5PushGroup(viewer, "/topology");
327:   ISView(orderIS, viewer);
328:   ISView(conesIS, viewer);
329:   ISView(cellsIS, viewer);
330:   ISView(orntsIS, viewer);
331:   PetscViewerHDF5PopGroup(viewer);
332:   ISDestroy(&orderIS);
333:   ISDestroy(&conesIS);
334:   ISDestroy(&cellsIS);
335:   ISDestroy(&orntsIS);
336:   ISRestoreIndices(globalPointNumbers, &gpoint);

338:   PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
339:   return(0);
340: }

342: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
343: {
344:   DMLabel         cutLabel;
345:   PetscSection    cSection;
346:   IS              cellIS, globalVertexNumbers;
347:   const PetscInt *gvertex;
348:   PetscInt       *vertices;
349:   IS              cutvertices;
350:   const PetscInt *cutverts;
351:   PetscInt        dim, depth, vStart, vEnd, vExtra = 0, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
352:   hid_t           fileId, groupId;
353:   PetscErrorCode  ierr;

356:   DMGetDimension(dm, &dim);
357:   DMPlexGetDepth(dm, &depth);
358:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
359:   DMGetCoordinateSection(dm, &cSection);
360:   DMPlexGetVTKCellHeight(dm, &cellHeight);
361:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
362:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
363:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
364:   for (cell = cStart; cell < cEnd; ++cell) {
365:     PetscInt *closure = NULL;
366:     PetscInt  closureSize, v, Nc = 0;

368:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
369:     for (v = 0; v < closureSize*2; v += 2) {
370:       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
371:     }
372:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
373:     conesSize += Nc;
374:     if (!numCornersLocal)           numCornersLocal = Nc;
375:     else if (numCornersLocal != Nc) numCornersLocal = 1;
376:   }
377:   MPIU_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
378:   if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");

380:   PetscViewerHDF5PushGroup(viewer, "/viz");
381:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
382:   PetscStackCallHDF5(H5Gclose,(groupId));
383:   PetscViewerHDF5PopGroup(viewer);

385:   DMGetLabel(dm, "periodic_cut", &cutLabel);
386:   if (cutLabel) {
387:     DMLabelGetStratumIS(cutLabel, 1, &cutvertices);
388:     ISGetIndices(cutvertices, &cutverts);
389:     ISGetLocalSize(cutvertices, &vExtra);
390:   }
391:   DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
392:   ISGetIndices(globalVertexNumbers, &gvertex);
393:   PetscMalloc1(conesSize, &vertices);
394:   for (cell = cStart, v = 0; cell < cEnd; ++cell) {
395:     PetscInt *closure = NULL;
396:     PetscInt  closureSize, Nc = 0, p;
397:     PetscBool replace = PETSC_FALSE;

399:     if (cutLabel) {
400:       PetscInt value;
401:       DMLabelGetValue(cutLabel, cell, &value);
402:       if (value == 2) replace = PETSC_TRUE;
403:     }
404:     DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
405:     for (p = 0; p < closureSize*2; p += 2) {
406:       if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
407:         closure[Nc++] = closure[p];
408:       }
409:     }
410:     DMPlexInvertCell_Internal(dim, Nc, closure);
411:     for (p = 0; p < Nc; ++p) {
412:       const PetscInt gv = gvertex[closure[p] - vStart];
413:       vertices[v++] = gv < 0 ? -(gv+1) : gv;
414:       if (replace) {
415:         PetscInt newv;
416:         PetscFindInt(closure[p], vExtra, cutverts, &newv);
417:         if (newv >= 0) vertices[v-1] = vEnd - vStart + newv;
418:       }
419:     }
420:     DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
421:   }
422:   if (cutLabel) {
423:     ISRestoreIndices(cutvertices, &cutverts);
424:     ISDestroy(&cutvertices);
425:   }
426:   if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
427:   ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);
428:   PetscLayoutSetBlockSize(cellIS->map, numCorners);
429:   PetscObjectSetName((PetscObject) cellIS, "cells");
430:   PetscViewerHDF5PushGroup(viewer, "/viz/topology");
431:   ISView(cellIS, viewer);
432:   PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
433:   ISDestroy(&cellIS);

435:   PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
436:   return(0);
437: }

439: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
440: {
441:   DM             cdm;
442:   Vec            coordinates, newcoords;
443:   PetscReal      lengthScale;
444:   PetscInt       m, M, bs;

448:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
449:   DMGetCoordinateDM(dm, &cdm);
450:   DMGetCoordinates(dm, &coordinates);
451:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
452:   PetscObjectSetName((PetscObject) newcoords, "vertices");
453:   VecGetSize(coordinates, &M);
454:   VecGetLocalSize(coordinates, &m);
455:   VecSetSizes(newcoords, m, M);
456:   VecGetBlockSize(coordinates, &bs);
457:   VecSetBlockSize(newcoords, bs);
458:   VecSetType(newcoords,VECSTANDARD);
459:   VecCopy(coordinates, newcoords);
460:   VecScale(newcoords, lengthScale);
461:   /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
462:   PetscViewerHDF5PushGroup(viewer, "/geometry");
463:   PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
464:   VecView(newcoords, viewer);
465:   PetscViewerPopFormat(viewer);
466:   PetscViewerHDF5PopGroup(viewer);
467:   VecDestroy(&newcoords);
468:   return(0);
469: }

471: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
472: {
473:   Vec              coordinates, newcoords;
474:   PetscSection     cSection;
475:   PetscScalar     *coords, *ncoords;
476:   DMLabel          cutLabel;
477:   const PetscReal *L;
478:   const DMBoundaryType *bd;
479:   PetscReal        lengthScale;
480:   PetscInt         vStart, vEnd, vExtra = 0, v, bs, coordSize, dof, off, d;
481:   PetscBool        localized, embedded;
482:   hid_t            fileId, groupId;
483:   PetscErrorCode   ierr;

486:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
487:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
488:   DMGetCoordinatesLocal(dm, &coordinates);
489:   VecGetBlockSize(coordinates, &bs);
490:   VecGetLocalSize(coordinates, &coordSize);
491:   DMGetCoordinatesLocalized(dm,&localized);
492:   if (localized == PETSC_FALSE) return(0);
493:   DMGetLabel(dm, "periodic_cut", &cutLabel);
494:   if (cutLabel) {
495:     IS              vertices;
496:     const PetscInt *verts;
497:     PetscInt        n;

499:     DMLabelGetStratumIS(cutLabel, 1, &vertices);
500:     ISGetIndices(vertices, &verts);
501:     ISGetLocalSize(vertices, &n);
502:     for (v = 0; v < n; ++v) {
503:       if ((verts[v] >= vStart) && (verts[v] < vEnd)) ++vExtra;
504:     }
505:     ISRestoreIndices(vertices, &verts);
506:     ISDestroy(&vertices);
507:   }
508:   DMGetPeriodicity(dm, NULL, NULL, &L, &bd);
509:   DMGetCoordinateSection(dm, &cSection);
510:   VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
511:   PetscSectionGetDof(cSection, vStart, &dof);
512:   embedded  = (PetscBool) (L && dof == 2 && !cutLabel);
513:   coordSize = 0;
514:   coordSize += dof*vExtra;
515:   for (v = vStart; v < vEnd; ++v) {
516:     PetscSectionGetDof(cSection, v, &dof);
517:     if (embedded) coordSize += dof+1;
518:     else          coordSize += dof;
519:   }
520:   if (embedded) {VecSetBlockSize(newcoords, bs+1);}
521:   else          {VecSetBlockSize(newcoords, bs);}
522:   VecSetSizes(newcoords, coordSize, PETSC_DETERMINE);
523:   VecSetType(newcoords, VECSTANDARD);
524:   VecGetArray(coordinates, &coords);
525:   VecGetArray(newcoords,   &ncoords);
526:   coordSize = 0;
527:   for (v = vStart; v < vEnd; ++v) {
528:     PetscSectionGetDof(cSection, v, &dof);
529:     PetscSectionGetOffset(cSection, v, &off);
530:     if (embedded) {
531:       if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {
532:         PetscReal theta, phi, r, R;
533:         /* XY-periodic */
534:         /* Suppose its an y-z circle, then
535:              \hat r = (0, cos(th), sin(th)) \hat x = (1, 0, 0)
536:            and the circle in that plane is
537:              \hat r cos(phi) + \hat x sin(phi) */
538:         theta = 2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1];
539:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
540:         r     = L[0]/(2.0*PETSC_PI * 2.0*L[1]);
541:         R     = L[1]/(2.0*PETSC_PI);
542:         ncoords[coordSize++] =  PetscSinReal(phi) * r;
543:         ncoords[coordSize++] = -PetscCosReal(theta) * (R + r * PetscCosReal(phi));
544:         ncoords[coordSize++] =  PetscSinReal(theta) * (R + r * PetscCosReal(phi));
545:       } else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
546:         /* X-periodic */
547:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
548:         ncoords[coordSize++] = coords[off+1];
549:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0])*(L[0]/(2.0*PETSC_PI));
550:       } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
551:         /* Y-periodic */
552:         ncoords[coordSize++] = coords[off+0];
553:         ncoords[coordSize++] = PetscSinReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
554:         ncoords[coordSize++] = -PetscCosReal(2.0*PETSC_PI*PetscRealPart(coords[off+1])/L[1])*(L[1]/(2.0*PETSC_PI));
555:       } else if ((bd[0] == DM_BOUNDARY_TWIST)) {
556:         PetscReal phi, r, R;
557:         /* Mobius strip */
558:         /* Suppose its an x-z circle, then
559:              \hat r = (-cos(phi), 0, sin(phi)) \hat y = (0, 1, 0)
560:            and in that plane we rotate by pi as we go around the circle
561:              \hat r cos(phi/2) + \hat y sin(phi/2) */
562:         phi   = 2.0*PETSC_PI*PetscRealPart(coords[off+0])/L[0];
563:         R     = L[0];
564:         r     = PetscRealPart(coords[off+1]) - L[1]/2.0;
565:         ncoords[coordSize++] = -PetscCosReal(phi) * (R + r * PetscCosReal(phi/2.0));
566:         ncoords[coordSize++] =  PetscSinReal(phi/2.0) * r;
567:         ncoords[coordSize++] =  PetscSinReal(phi) * (R + r * PetscCosReal(phi/2.0));
568:       } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
569:     } else {
570:       if (cutLabel) {
571:         DMLocalizeCoordinate(dm, &coords[off], PETSC_TRUE, &ncoords[coordSize]);
572:         coordSize += dof;
573:       } else {
574:         for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
575:       }
576:     }
577:   }
578:   if (cutLabel) {
579:     IS              vertices;
580:     const PetscInt *verts;
581:     PetscInt        n;

583:     DMLabelGetStratumIS(cutLabel, 1, &vertices);
584:     ISGetIndices(vertices, &verts);
585:     ISGetLocalSize(vertices, &n);
586:     for (v = 0; v < n; ++v) {
587:       if ((verts[v] < vStart) || (verts[v] >= vEnd)) continue;
588:       PetscSectionGetDof(cSection, verts[v], &dof);
589:       PetscSectionGetOffset(cSection, verts[v], &off);
590:       for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d] + ((bd[d] == DM_BOUNDARY_PERIODIC) ? L[d] : 0.0);
591:     }
592:     ISRestoreIndices(vertices, &verts);
593:     ISDestroy(&vertices);
594:   }
595:   VecRestoreArray(coordinates, &coords);
596:   VecRestoreArray(newcoords,   &ncoords);
597:   PetscObjectSetName((PetscObject) newcoords, "vertices");
598:   VecScale(newcoords, lengthScale);
599:   PetscViewerHDF5PushGroup(viewer, "/viz");
600:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
601:   PetscStackCallHDF5(H5Gclose,(groupId));
602:   PetscViewerHDF5PopGroup(viewer);
603:   PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
604:   VecView(newcoords, viewer);
605:   PetscViewerHDF5PopGroup(viewer);
606:   VecDestroy(&newcoords);
607:   return(0);
608: }


611: static PetscErrorCode DMPlexWriteLabels_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
612: {
613:   const PetscInt   *gpoint;
614:   PetscInt          numLabels, l;
615:   hid_t             fileId, groupId;
616:   PetscErrorCode    ierr;

619:   ISGetIndices(globalPointNumbers, &gpoint);
620:   PetscViewerHDF5PushGroup(viewer, "/labels");
621:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
622:   if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
623:   PetscViewerHDF5PopGroup(viewer);
624:   DMGetNumLabels(dm, &numLabels);
625:   for (l = 0; l < numLabels; ++l) {
626:     DMLabel         label;
627:     const char     *name;
628:     IS              valueIS, pvalueIS, globalValueIS;
629:     const PetscInt *values;
630:     PetscInt        numValues, v;
631:     PetscBool       isDepth, output;
632:     char            group[PETSC_MAX_PATH_LEN];

634:     DMGetLabelName(dm, l, &name);
635:     DMGetLabelOutput(dm, name, &output);
636:     PetscStrncmp(name, "depth", 10, &isDepth);
637:     if (isDepth || !output) continue;
638:     DMGetLabel(dm, name, &label);
639:     PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
640:     PetscViewerHDF5PushGroup(viewer, group);
641:     PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
642:     if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
643:     PetscViewerHDF5PopGroup(viewer);
644:     DMLabelGetValueIS(label, &valueIS);
645:     /* Must copy to a new IS on the global comm */
646:     ISGetLocalSize(valueIS, &numValues);
647:     ISGetIndices(valueIS, &values);
648:     ISCreateGeneral(PetscObjectComm((PetscObject) dm), numValues, values, PETSC_COPY_VALUES, &pvalueIS);
649:     ISRestoreIndices(valueIS, &values);
650:     ISAllGather(pvalueIS, &globalValueIS);
651:     ISDestroy(&pvalueIS);
652:     ISSortRemoveDups(globalValueIS);
653:     ISGetLocalSize(globalValueIS, &numValues);
654:     ISGetIndices(globalValueIS, &values);
655:     for (v = 0; v < numValues; ++v) {
656:       IS              stratumIS, globalStratumIS;
657:       const PetscInt *spoints = NULL;
658:       PetscInt       *gspoints, n = 0, gn, p;
659:       const char     *iname = "indices";

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

664:       if (stratumIS) {ISGetLocalSize(stratumIS, &n);}
665:       if (stratumIS) {ISGetIndices(stratumIS, &spoints);}
666:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
667:       PetscMalloc1(gn,&gspoints);
668:       for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
669:       if (stratumIS) {ISRestoreIndices(stratumIS, &spoints);}
670:       ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
671:       if (stratumIS) {PetscObjectGetName((PetscObject) stratumIS, &iname);}
672:       PetscObjectSetName((PetscObject) globalStratumIS, iname);

674:       PetscViewerHDF5PushGroup(viewer, group);
675:       ISView(globalStratumIS, viewer);
676:       PetscViewerHDF5PopGroup(viewer);
677:       ISDestroy(&globalStratumIS);
678:       ISDestroy(&stratumIS);
679:     }
680:     ISRestoreIndices(globalValueIS, &values);
681:     ISDestroy(&globalValueIS);
682:     ISDestroy(&valueIS);
683:   }
684:   ISRestoreIndices(globalPointNumbers, &gpoint);
685:   return(0);
686: }

688: /* We only write cells and vertices. Does this screw up parallel reading? */
689: PetscErrorCode DMPlexView_HDF5_Internal(DM dm, PetscViewer viewer)
690: {
691:   IS                globalPointNumbers;
692:   PetscViewerFormat format;
693:   PetscBool         viz_geom=PETSC_FALSE, xdmf_topo=PETSC_FALSE, petsc_topo=PETSC_FALSE;
694:   PetscErrorCode    ierr;

697:   DMPlexCreatePointNumbering(dm, &globalPointNumbers);
698:   DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
699:   DMPlexWriteLabels_HDF5_Static(dm, globalPointNumbers, viewer);

701:   PetscViewerGetFormat(viewer, &format);
702:   switch (format) {
703:     case PETSC_VIEWER_HDF5_VIZ:
704:       viz_geom    = PETSC_TRUE;
705:       xdmf_topo   = PETSC_TRUE;
706:       break;
707:     case PETSC_VIEWER_HDF5_XDMF:
708:       xdmf_topo   = PETSC_TRUE;
709:       break;
710:     case PETSC_VIEWER_HDF5_PETSC:
711:       petsc_topo  = PETSC_TRUE;
712:       break;
713:     case PETSC_VIEWER_DEFAULT:
714:     case PETSC_VIEWER_NATIVE:
715:       viz_geom    = PETSC_TRUE;
716:       xdmf_topo   = PETSC_TRUE;
717:       petsc_topo  = PETSC_TRUE;
718:       break;
719:     default:
720:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "PetscViewerFormat %s not supported for HDF5 output.", PetscViewerFormats[format]);
721:   }

723:   if (viz_geom)   {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
724:   if (xdmf_topo)  {DMPlexWriteTopology_Vertices_HDF5_Static(dm, viewer);}
725:   if (petsc_topo) {DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);}

727:   ISDestroy(&globalPointNumbers);
728:   return(0);
729: }

731: typedef struct {
732:   PetscMPIInt rank;
733:   DM          dm;
734:   PetscViewer viewer;
735:   DMLabel     label;
736: } LabelCtx;

738: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
739: {
740:   PetscViewer     viewer = ((LabelCtx *) op_data)->viewer;
741:   DMLabel         label  = ((LabelCtx *) op_data)->label;
742:   IS              stratumIS;
743:   const PetscInt *ind;
744:   PetscInt        value, N, i;
745:   const char     *lname;
746:   char            group[PETSC_MAX_PATH_LEN];
747:   PetscErrorCode  ierr;

749:   PetscOptionsStringToInt(name, &value);
750:   ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
751:   PetscObjectSetName((PetscObject) stratumIS, "indices");
752:   DMLabelGetName(label, &lname);
753:   PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
754:   PetscViewerHDF5PushGroup(viewer, group);
755:   {
756:     /* Force serial load */
757:     PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
758:     PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
759:     PetscLayoutSetSize(stratumIS->map, N);
760:   }
761:   ISLoad(stratumIS, viewer);
762:   PetscViewerHDF5PopGroup(viewer);
763:   ISGetLocalSize(stratumIS, &N);
764:   ISGetIndices(stratumIS, &ind);
765:   for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
766:   ISRestoreIndices(stratumIS, &ind);
767:   ISDestroy(&stratumIS);
768:   return 0;
769: }

771: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
772: {
773:   DM             dm  = ((LabelCtx *) op_data)->dm;
774:   hsize_t        idx = 0;
776:   herr_t         err;

778:   DMCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
779:   DMGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
780:   PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
781:   return err;
782: }

784: PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM dm, PetscViewer viewer)
785: {
786:   LabelCtx        ctx;
787:   hid_t           fileId, groupId;
788:   hsize_t         idx = 0;
789:   PetscErrorCode  ierr;

792:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &ctx.rank);
793:   ctx.dm     = dm;
794:   ctx.viewer = viewer;
795:   PetscViewerHDF5PushGroup(viewer, "/labels");
796:   PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
797:   PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
798:   PetscStackCallHDF5(H5Gclose,(groupId));
799:   PetscViewerHDF5PopGroup(viewer);
800:   return(0);
801: }

803: /* The first version will read everything onto proc 0, letting the user distribute
804:    The next will create a naive partition, and then rebalance after reading
805: */
806: PetscErrorCode DMPlexLoad_HDF5_Internal(DM dm, PetscViewer viewer)
807: {
808:   PetscSection    coordSection;
809:   Vec             coordinates;
810:   IS              orderIS, conesIS, cellsIS, orntsIS;
811:   const PetscInt *order, *cones, *cells, *ornts;
812:   PetscReal       lengthScale;
813:   PetscInt       *cone, *ornt;
814:   PetscInt        dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
815:   PetscMPIInt     rank;
816:   PetscErrorCode  ierr;

819:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
820:   /* Read toplogy */
821:   PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
822:   DMSetDimension(dm, dim);
823:   PetscViewerHDF5PushGroup(viewer, "/topology");

825:   ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
826:   PetscObjectSetName((PetscObject) orderIS, "order");
827:   ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
828:   PetscObjectSetName((PetscObject) conesIS, "cones");
829:   ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
830:   PetscObjectSetName((PetscObject) cellsIS, "cells");
831:   ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
832:   PetscObjectSetName((PetscObject) orntsIS, "orientation");
833:   {
834:     /* Force serial load */
835:     PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
836:     PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
837:     PetscLayoutSetSize(orderIS->map, pEnd);
838:     PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
839:     PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
840:     PetscLayoutSetSize(conesIS->map, pEnd);
841:     pEnd = !rank ? pEnd : 0;
842:     PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
843:     PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
844:     PetscLayoutSetSize(cellsIS->map, N);
845:     PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
846:     PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
847:     PetscLayoutSetSize(orntsIS->map, N);
848:   }
849:   ISLoad(orderIS, viewer);
850:   ISLoad(conesIS, viewer);
851:   ISLoad(cellsIS, viewer);
852:   ISLoad(orntsIS, viewer);
853:   PetscViewerHDF5PopGroup(viewer);
854:   /* Read geometry */
855:   PetscViewerHDF5PushGroup(viewer, "/geometry");
856:   VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
857:   PetscObjectSetName((PetscObject) coordinates, "vertices");
858:   {
859:     /* Force serial load */
860:     PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
861:     VecSetSizes(coordinates, !rank ? N : 0, N);
862:     VecSetBlockSize(coordinates, spatialDim);
863:   }
864:   VecLoad(coordinates, viewer);
865:   PetscViewerHDF5PopGroup(viewer);
866:   DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
867:   VecScale(coordinates, 1.0/lengthScale);
868:   VecGetLocalSize(coordinates, &numVertices);
869:   VecGetBlockSize(coordinates, &spatialDim);
870:   numVertices /= spatialDim;
871:   /* Create Plex */
872:   DMPlexSetChart(dm, 0, pEnd);
873:   ISGetIndices(orderIS, &order);
874:   ISGetIndices(conesIS, &cones);
875:   for (p = 0; p < pEnd; ++p) {
876:     DMPlexSetConeSize(dm, order[p], cones[p]);
877:     maxConeSize = PetscMax(maxConeSize, cones[p]);
878:   }
879:   DMSetUp(dm);
880:   ISGetIndices(cellsIS, &cells);
881:   ISGetIndices(orntsIS, &ornts);
882:   PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
883:   for (p = 0, q = 0; p < pEnd; ++p) {
884:     for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
885:     DMPlexSetCone(dm, order[p], cone);
886:     DMPlexSetConeOrientation(dm, order[p], ornt);
887:   }
888:   PetscFree2(cone,ornt);
889:   ISRestoreIndices(orderIS, &order);
890:   ISRestoreIndices(conesIS, &cones);
891:   ISRestoreIndices(cellsIS, &cells);
892:   ISRestoreIndices(orntsIS, &ornts);
893:   ISDestroy(&orderIS);
894:   ISDestroy(&conesIS);
895:   ISDestroy(&cellsIS);
896:   ISDestroy(&orntsIS);
897:   DMPlexSymmetrize(dm);
898:   DMPlexStratify(dm);
899:   /* Create coordinates */
900:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
901:   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);
902:   DMGetCoordinateSection(dm, &coordSection);
903:   PetscSectionSetNumFields(coordSection, 1);
904:   PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
905:   PetscSectionSetChart(coordSection, vStart, vEnd);
906:   for (v = vStart; v < vEnd; ++v) {
907:     PetscSectionSetDof(coordSection, v, spatialDim);
908:     PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
909:   }
910:   PetscSectionSetUp(coordSection);
911:   DMSetCoordinates(dm, coordinates);
912:   VecDestroy(&coordinates);
913:   /* Read Labels */
914:   DMPlexLoadLabels_HDF5_Internal(dm, viewer);
915:   return(0);
916: }
917: #endif