Actual source code: plexhdf5.c
petsc-3.7.7 2017-09-25
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, §ion);
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, §ionGlobal);
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