Actual source code: plexhdf5.c
petsc-3.6.1 2015-08-06
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_Seq(Vec, PetscViewer);
7: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
9: #if defined(PETSC_HAVE_HDF5)
12: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
13: {
14: Vec stamp;
15: PetscMPIInt rank;
19: if (seqnum < 0) return(0);
20: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
21: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
22: VecSetBlockSize(stamp, 1);
23: PetscObjectSetName((PetscObject) stamp, seqname);
24: if (!rank) {
25: PetscReal timeScale;
26: PetscBool istime;
28: PetscStrncmp(seqname, "time", 5, &istime);
29: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); value *= timeScale;}
30: VecSetValue(stamp, 0, value, INSERT_VALUES);
31: }
32: VecAssemblyBegin(stamp);
33: VecAssemblyEnd(stamp);
34: PetscViewerHDF5PushGroup(viewer, "/");
35: PetscViewerHDF5SetTimestep(viewer, seqnum);
36: VecView(stamp, viewer);
37: PetscViewerHDF5PopGroup(viewer);
38: VecDestroy(&stamp);
39: return(0);
40: }
44: PetscErrorCode DMSequenceLoad_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
45: {
46: Vec stamp;
47: PetscMPIInt rank;
51: if (seqnum < 0) return(0);
52: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
53: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
54: VecSetBlockSize(stamp, 1);
55: PetscObjectSetName((PetscObject) stamp, seqname);
56: PetscViewerHDF5PushGroup(viewer, "/");
57: PetscViewerHDF5SetTimestep(viewer, seqnum);
58: VecLoad(stamp, viewer);
59: PetscViewerHDF5PopGroup(viewer);
60: if (!rank) {
61: const PetscScalar *a;
62: PetscReal timeScale;
63: PetscBool istime;
65: VecGetArrayRead(stamp, &a);
66: *value = a[0];
67: VecRestoreArrayRead(stamp, &a);
68: PetscStrncmp(seqname, "time", 5, &istime);
69: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
70: }
71: VecDestroy(&stamp);
72: return(0);
73: }
77: PetscErrorCode VecView_Plex_Local_HDF5(Vec v, PetscViewer viewer)
78: {
79: DM dm;
80: DM dmBC;
81: PetscSection section, sectionGlobal;
82: Vec gv;
83: const char *name;
84: PetscViewerVTKFieldType ft;
85: PetscViewerFormat format;
86: PetscInt seqnum;
87: PetscReal seqval;
88: PetscBool isseq;
89: PetscErrorCode ierr;
92: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
93: VecGetDM(v, &dm);
94: DMGetDefaultSection(dm, §ion);
95: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
96: PetscViewerHDF5SetTimestep(viewer, seqnum);
97: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
98: PetscViewerGetFormat(viewer, &format);
99: DMGetOutputDM(dm, &dmBC);
100: DMGetDefaultGlobalSection(dmBC, §ionGlobal);
101: DMGetGlobalVector(dmBC, &gv);
102: PetscObjectGetName((PetscObject) v, &name);
103: PetscObjectSetName((PetscObject) gv, name);
104: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
105: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
106: PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
107: if (isseq) {VecView_Seq(gv, viewer);}
108: else {VecView_MPI(gv, viewer);}
109: if (format == PETSC_VIEWER_HDF5_VIZ) {
110: /* Output visualization representation */
111: PetscInt numFields, f;
113: PetscSectionGetNumFields(section, &numFields);
114: for (f = 0; f < numFields; ++f) {
115: Vec subv;
116: IS is;
117: const char *fname, *fgroup;
118: char subname[PETSC_MAX_PATH_LEN];
119: char group[PETSC_MAX_PATH_LEN];
120: PetscInt pStart, pEnd;
121: PetscBool flag;
123: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
124: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
125: PetscSectionGetFieldName(section, f, &fname);
126: if (!fname) continue;
127: PetscViewerHDF5PushGroup(viewer, fgroup);
128: PetscSectionGetField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
129: PetscStrcpy(subname, name);
130: PetscStrcat(subname, "_");
131: PetscStrcat(subname, fname);
132: PetscObjectSetName((PetscObject) subv, subname);
133: if (isseq) {VecView_Seq(subv, viewer);}
134: else {VecView_MPI(subv, viewer);}
135: PetscSectionRestoreField_Internal(section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
136: PetscViewerHDF5PopGroup(viewer);
137: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s/%s", fgroup, subname);
138: PetscViewerHDF5HasAttribute(viewer, group, "vector_field_type", &flag);
139: if (!flag) {
140: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
141: PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "vector");
142: } else {
143: PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "scalar");
144: }
145: }
146: }
147: }
148: DMRestoreGlobalVector(dmBC, &gv);
149: return(0);
150: }
154: PetscErrorCode VecView_Plex_HDF5(Vec v, PetscViewer viewer)
155: {
156: DM dm;
157: Vec locv;
158: const char *name;
162: VecGetDM(v, &dm);
163: DMGetLocalVector(dm, &locv);
164: PetscObjectGetName((PetscObject) v, &name);
165: PetscObjectSetName((PetscObject) locv, name);
166: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
167: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
168: PetscViewerHDF5PushGroup(viewer, "/fields");
169: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
170: VecView_Plex_Local(locv, viewer);
171: PetscViewerPopFormat(viewer);
172: PetscViewerHDF5PopGroup(viewer);
173: DMRestoreLocalVector(dm, &locv);
174: return(0);
175: }
179: PetscErrorCode VecLoad_Plex_HDF5(Vec v, PetscViewer viewer)
180: {
181: DM dm;
182: Vec locv;
183: const char *name;
184: PetscInt seqnum;
188: VecGetDM(v, &dm);
189: DMGetLocalVector(dm, &locv);
190: PetscObjectGetName((PetscObject) v, &name);
191: PetscObjectSetName((PetscObject) locv, name);
192: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
193: PetscViewerHDF5SetTimestep(viewer, seqnum);
194: PetscViewerHDF5PushGroup(viewer, "/fields");
195: VecLoad_Plex_Local(locv, viewer);
196: PetscViewerHDF5PopGroup(viewer);
197: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
198: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
199: DMRestoreLocalVector(dm, &locv);
200: return(0);
201: }
205: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
206: {
207: IS orderIS, conesIS, cellsIS, orntsIS;
208: const PetscInt *gpoint;
209: PetscInt *order, *sizes, *cones, *ornts;
210: PetscInt dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
211: PetscErrorCode ierr;
214: ISGetIndices(globalPointNumbers, &gpoint);
215: DMGetDimension(dm, &dim);
216: DMPlexGetChart(dm, &pStart, &pEnd);
217: for (p = pStart; p < pEnd; ++p) {
218: if (gpoint[p] >= 0) {
219: PetscInt coneSize;
221: DMPlexGetConeSize(dm, p, &coneSize);
222: conesSize += 1;
223: cellsSize += coneSize;
224: }
225: }
226: PetscMalloc1(conesSize, &order);
227: PetscMalloc1(conesSize, &sizes);
228: PetscMalloc1(cellsSize, &cones);
229: PetscMalloc1(cellsSize, &ornts);
230: for (p = pStart; p < pEnd; ++p) {
231: if (gpoint[p] >= 0) {
232: const PetscInt *cone, *ornt;
233: PetscInt coneSize, cp;
235: DMPlexGetConeSize(dm, p, &coneSize);
236: DMPlexGetCone(dm, p, &cone);
237: DMPlexGetConeOrientation(dm, p, &ornt);
238: order[s] = gpoint[p];
239: sizes[s++] = coneSize;
240: for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
241: }
242: }
243: if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %d != %d", s, conesSize);
244: if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %d != %d", c, cellsSize);
245: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
246: PetscObjectSetName((PetscObject) orderIS, "order");
247: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
248: PetscObjectSetName((PetscObject) conesIS, "cones");
249: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
250: PetscObjectSetName((PetscObject) cellsIS, "cells");
251: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
252: PetscObjectSetName((PetscObject) orntsIS, "orientation");
253: PetscViewerHDF5PushGroup(viewer, "/topology");
254: ISView(orderIS, viewer);
255: ISView(conesIS, viewer);
256: ISView(cellsIS, viewer);
257: ISView(orntsIS, viewer);
258: PetscViewerHDF5PopGroup(viewer);
259: ISDestroy(&orderIS);
260: ISDestroy(&conesIS);
261: ISDestroy(&cellsIS);
262: ISDestroy(&orntsIS);
263: ISRestoreIndices(globalPointNumbers, &gpoint);
265: PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
266: return(0);
267: }
271: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, DMLabel label, PetscInt labelId, PetscViewer viewer)
272: {
273: IS cellIS, globalVertexNumbers;
274: const PetscInt *gvertex;
275: PetscInt *vertices;
276: PetscInt dim, depth, vStart, vEnd, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
277: hid_t fileId, groupId;
278: PetscErrorCode ierr;
281: DMGetDimension(dm, &dim);
282: DMPlexGetDepth(dm, &depth);
283: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
284: DMPlexGetVTKCellHeight(dm, &cellHeight);
285: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
286: DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, PETSC_NULL);
287: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
288: if (depth == 1) return(0);
289: for (cell = cStart; cell < cEnd; ++cell) {
290: PetscInt *closure = NULL;
291: PetscInt closureSize, v, Nc = 0;
293: if (label) {
294: PetscInt value;
295: DMLabelGetValue(label, cell, &value);
296: if (value == labelId) continue;
297: }
298: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
299: for (v = 0; v < closureSize*2; v += 2) {
300: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
301: }
302: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
303: conesSize += Nc;
304: if (!numCornersLocal) numCornersLocal = Nc;
305: else if (numCornersLocal != Nc) numCornersLocal = 1;
306: }
307: MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
308: if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
310: PetscViewerHDF5PushGroup(viewer, "/viz");
311: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
312: PetscStackCallHDF5(H5Gclose,(groupId));
313: PetscViewerHDF5PopGroup(viewer);
315: DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
316: ISGetIndices(globalVertexNumbers, &gvertex);
317: PetscMalloc1(conesSize, &vertices);
318: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
319: PetscInt *closure = NULL;
320: PetscInt closureSize, Nc = 0, p;
322: if (label) {
323: PetscInt value;
324: DMLabelGetValue(label, cell, &value);
325: if (value == labelId) continue;
326: }
327: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
328: for (p = 0; p < closureSize*2; p += 2) {
329: if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
330: closure[Nc++] = closure[p];
331: }
332: }
333: DMPlexInvertCell_Internal(dim, Nc, closure);
334: for (p = 0; p < Nc; ++p) {
335: const PetscInt gv = gvertex[closure[p] - vStart];
336: vertices[v++] = gv < 0 ? -(gv+1) : gv;
337: }
338: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
339: }
340: if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
341: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);
342: PetscLayoutSetBlockSize(cellIS->map, numCorners);
343: PetscObjectSetName((PetscObject) cellIS, "cells");
344: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
345: ISView(cellIS, viewer);
346: #if 0
347: if (numCorners == 1) {
348: ISView(coneIS, viewer);
349: } else {
350: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
351: }
352: ISDestroy(&coneIS);
353: #else
354: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
355: #endif
356: ISDestroy(&cellIS);
358: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
359: return(0);
360: }
364: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
365: {
366: DM cdm;
367: Vec coordinates, newcoords;
368: PetscReal lengthScale;
369: PetscInt m, M, bs;
373: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
374: DMGetCoordinateDM(dm, &cdm);
375: DMGetCoordinates(dm, &coordinates);
376: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
377: PetscObjectSetName((PetscObject) newcoords, "vertices");
378: VecGetSize(coordinates, &M);
379: VecGetLocalSize(coordinates, &m);
380: VecSetSizes(newcoords, m, M);
381: VecGetBlockSize(coordinates, &bs);
382: VecSetBlockSize(newcoords, bs);
383: VecSetType(newcoords,VECSTANDARD);
384: VecCopy(coordinates, newcoords);
385: VecScale(newcoords, lengthScale);
386: /* Did not use DMGetGlobalVector() in order to bypass default group assignment */
387: PetscViewerHDF5PushGroup(viewer, "/geometry");
388: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
389: VecView(newcoords, viewer);
390: PetscViewerPopFormat(viewer);
391: PetscViewerHDF5PopGroup(viewer);
392: VecDestroy(&newcoords);
393: return(0);
394: }
398: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
399: {
400: Vec coordinates, newcoords;
401: PetscSection cSection;
402: PetscScalar *coords, *ncoords;
403: const PetscReal *L;
404: const DMBoundaryType *bd;
405: PetscReal lengthScale;
406: PetscInt vStart, vEnd, v, bs, coordSize, dof, off, d;
407: hid_t fileId, groupId;
408: PetscErrorCode ierr;
411: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
412: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
413: DMGetCoordinatesLocal(dm, &coordinates);
414: VecGetBlockSize(coordinates, &bs);
415: VecGetLocalSize(coordinates, &coordSize);
416: if (coordSize == (vEnd - vStart)*bs) return(0);
417: DMGetPeriodicity(dm, NULL, &L, &bd);
418: DMGetCoordinateSection(dm, &cSection);
419: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
420: coordSize = 0;
421: for (v = vStart; v < vEnd; ++v) {
422: PetscSectionGetDof(cSection, v, &dof);
423: if (L && dof == 2) coordSize += dof+1;
424: else coordSize += dof;
425: }
426: if (L && bs == 2) {VecSetBlockSize(newcoords, bs+1);}
427: else {VecSetBlockSize(newcoords, bs);}
428: VecSetSizes(newcoords, coordSize, PETSC_DETERMINE);
429: VecSetType(newcoords,VECSTANDARD);
430: VecGetArray(coordinates, &coords);
431: VecGetArray(newcoords, &ncoords);
432: coordSize = 0;
433: for (v = vStart; v < vEnd; ++v) {
434: PetscSectionGetDof(cSection, v, &dof);
435: PetscSectionGetOffset(cSection, v, &off);
436: if (L && dof == 2) {
437: /* Need to do torus */
438: if ((bd[0] == DM_BOUNDARY_PERIODIC) && (bd[1] == DM_BOUNDARY_PERIODIC)) {SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot embed doubly periodic domain yet");}
439: else if ((bd[0] == DM_BOUNDARY_PERIODIC)) {
440: /* X-periodic */
441: ncoords[coordSize++] = -cos(2.0*PETSC_PI*coords[off+0]/L[0])*(L[0]/(2.0*PETSC_PI));
442: ncoords[coordSize++] = coords[off+1];
443: ncoords[coordSize++] = sin(2.0*PETSC_PI*coords[off+0]/L[0])*(L[0]/(2.0*PETSC_PI));
444: } else if ((bd[1] == DM_BOUNDARY_PERIODIC)) {
445: /* Y-periodic */
446: ncoords[coordSize++] = coords[off+0];
447: ncoords[coordSize++] = sin(2.0*PETSC_PI*coords[off+1]/L[1])*(L[1]/(2.0*PETSC_PI));
448: ncoords[coordSize++] = -cos(2.0*PETSC_PI*coords[off+1]/L[1])*(L[1]/(2.0*PETSC_PI));
449: } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Cannot handle periodicity in this domain");
450: } else {
451: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
452: }
453: }
454: VecRestoreArray(coordinates, &coords);
455: VecRestoreArray(newcoords, &ncoords);
456: PetscObjectSetName((PetscObject) newcoords, "vertices");
457: VecScale(newcoords, lengthScale);
458: PetscViewerHDF5PushGroup(viewer, "/viz");
459: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
460: PetscStackCallHDF5(H5Gclose,(groupId));
461: PetscViewerHDF5PopGroup(viewer);
462: PetscViewerHDF5PushGroup(viewer, "/viz/geometry");
463: VecView(newcoords, viewer);
464: PetscViewerHDF5PopGroup(viewer);
465: VecDestroy(&newcoords);
466: return(0);
467: }
471: /* We only write cells and vertices. Does this screw up parallel reading? */
472: PetscErrorCode DMPlexView_HDF5(DM dm, PetscViewer viewer)
473: {
474: DMLabel label = NULL;
475: PetscInt labelId = 0;
476: IS globalPointNumbers;
477: const PetscInt *gpoint;
478: PetscInt numLabels, l;
479: hid_t fileId, groupId;
480: PetscViewerFormat format;
481: PetscErrorCode ierr;
484: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
485: PetscViewerGetFormat(viewer, &format);
486: DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
487: if (format == PETSC_VIEWER_HDF5_VIZ) {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
488: DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);
489: if (format == PETSC_VIEWER_HDF5_VIZ) {DMPlexWriteTopology_Vertices_HDF5_Static(dm, label, labelId, viewer);}
490: /* Write Labels*/
491: ISGetIndices(globalPointNumbers, &gpoint);
492: PetscViewerHDF5PushGroup(viewer, "/labels");
493: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
494: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
495: PetscViewerHDF5PopGroup(viewer);
496: DMPlexGetNumLabels(dm, &numLabels);
497: for (l = 0; l < numLabels; ++l) {
498: DMLabel label;
499: const char *name;
500: IS valueIS, globalValueIS;
501: const PetscInt *values;
502: PetscInt numValues, v;
503: PetscBool isDepth, output;
504: char group[PETSC_MAX_PATH_LEN];
506: DMPlexGetLabelName(dm, l, &name);
507: DMPlexGetLabelOutput(dm, name, &output);
508: PetscStrncmp(name, "depth", 10, &isDepth);
509: if (isDepth || !output) continue;
510: DMPlexGetLabel(dm, name, &label);
511: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s", name);
512: PetscViewerHDF5PushGroup(viewer, group);
513: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
514: if (groupId != fileId) PetscStackCallHDF5(H5Gclose,(groupId));
515: PetscViewerHDF5PopGroup(viewer);
516: DMLabelGetValueIS(label, &valueIS);
517: ISAllGather(valueIS, &globalValueIS);
518: ISSortRemoveDups(globalValueIS);
519: ISGetLocalSize(globalValueIS, &numValues);
520: ISGetIndices(globalValueIS, &values);
521: for (v = 0; v < numValues; ++v) {
522: IS stratumIS, globalStratumIS;
523: const PetscInt *spoints;
524: PetscInt *gspoints, n, gn, p;
525: const char *iname;
527: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%d", name, values[v]);
528: DMLabelGetStratumIS(label, values[v], &stratumIS);
530: ISGetLocalSize(stratumIS, &n);
531: ISGetIndices(stratumIS, &spoints);
532: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
533: PetscMalloc1(gn,&gspoints);
534: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
535: ISRestoreIndices(stratumIS, &spoints);
536: ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
537: PetscObjectGetName((PetscObject) stratumIS, &iname);
538: PetscObjectSetName((PetscObject) globalStratumIS, iname);
540: PetscViewerHDF5PushGroup(viewer, group);
541: ISView(globalStratumIS, viewer);
542: PetscViewerHDF5PopGroup(viewer);
543: ISDestroy(&globalStratumIS);
544: ISDestroy(&stratumIS);
545: }
546: ISRestoreIndices(globalValueIS, &values);
547: ISDestroy(&globalValueIS);
548: ISDestroy(&valueIS);
549: }
550: ISRestoreIndices(globalPointNumbers, &gpoint);
551: ISDestroy(&globalPointNumbers);
552: return(0);
553: }
555: typedef struct {
556: PetscMPIInt rank;
557: DM dm;
558: PetscViewer viewer;
559: DMLabel label;
560: } LabelCtx;
562: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
563: {
564: PetscViewer viewer = ((LabelCtx *) op_data)->viewer;
565: DMLabel label = ((LabelCtx *) op_data)->label;
566: IS stratumIS;
567: const PetscInt *ind;
568: PetscInt value, N, i;
569: const char *lname;
570: char group[PETSC_MAX_PATH_LEN];
571: PetscErrorCode ierr;
573: PetscOptionsStringToInt(name, &value);
574: ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
575: PetscObjectSetName((PetscObject) stratumIS, "indices");
576: DMLabelGetName(label, &lname);
577: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
578: PetscViewerHDF5PushGroup(viewer, group);
579: {
580: /* Force serial load */
581: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
582: PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
583: PetscLayoutSetSize(stratumIS->map, N);
584: }
585: ISLoad(stratumIS, viewer);
586: PetscViewerHDF5PopGroup(viewer);
587: ISGetLocalSize(stratumIS, &N);
588: ISGetIndices(stratumIS, &ind);
589: for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
590: ISRestoreIndices(stratumIS, &ind);
591: ISDestroy(&stratumIS);
592: return 0;
593: }
595: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
596: {
597: DM dm = ((LabelCtx *) op_data)->dm;
598: hsize_t idx = 0;
600: herr_t err;
602: DMPlexCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
603: DMPlexGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
604: PetscStackCall("H5Literate_by_name",err = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0));
605: return err;
606: }
610: /* The first version will read everything onto proc 0, letting the user distribute
611: The next will create a naive partition, and then rebalance after reading
612: */
613: PetscErrorCode DMPlexLoad_HDF5(DM dm, PetscViewer viewer)
614: {
615: LabelCtx ctx;
616: PetscSection coordSection;
617: Vec coordinates;
618: IS orderIS, conesIS, cellsIS, orntsIS;
619: const PetscInt *order, *cones, *cells, *ornts;
620: PetscReal lengthScale;
621: PetscInt *cone, *ornt;
622: PetscInt dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
623: hid_t fileId, groupId;
624: hsize_t idx = 0;
625: PetscMPIInt rank;
626: PetscErrorCode ierr;
629: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
630: /* Read toplogy */
631: PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
632: DMSetDimension(dm, dim);
633: PetscViewerHDF5PushGroup(viewer, "/topology");
635: ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
636: PetscObjectSetName((PetscObject) orderIS, "order");
637: ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
638: PetscObjectSetName((PetscObject) conesIS, "cones");
639: ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
640: PetscObjectSetName((PetscObject) cellsIS, "cells");
641: ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
642: PetscObjectSetName((PetscObject) orntsIS, "orientation");
643: {
644: /* Force serial load */
645: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
646: PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
647: PetscLayoutSetSize(orderIS->map, pEnd);
648: PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
649: PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
650: PetscLayoutSetSize(conesIS->map, pEnd);
651: pEnd = !rank ? pEnd : 0;
652: PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
653: PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
654: PetscLayoutSetSize(cellsIS->map, N);
655: PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
656: PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
657: PetscLayoutSetSize(orntsIS->map, N);
658: }
659: ISLoad(orderIS, viewer);
660: ISLoad(conesIS, viewer);
661: ISLoad(cellsIS, viewer);
662: ISLoad(orntsIS, viewer);
663: PetscViewerHDF5PopGroup(viewer);
664: /* Read geometry */
665: PetscViewerHDF5PushGroup(viewer, "/geometry");
666: VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
667: PetscObjectSetName((PetscObject) coordinates, "vertices");
668: {
669: /* Force serial load */
670: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
671: VecSetSizes(coordinates, !rank ? N : 0, N);
672: VecSetBlockSize(coordinates, spatialDim);
673: }
674: VecLoad(coordinates, viewer);
675: PetscViewerHDF5PopGroup(viewer);
676: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
677: VecScale(coordinates, 1.0/lengthScale);
678: VecGetLocalSize(coordinates, &numVertices);
679: VecGetBlockSize(coordinates, &spatialDim);
680: numVertices /= spatialDim;
681: /* Create Plex */
682: DMPlexSetChart(dm, 0, pEnd);
683: ISGetIndices(orderIS, &order);
684: ISGetIndices(conesIS, &cones);
685: for (p = 0; p < pEnd; ++p) {
686: DMPlexSetConeSize(dm, order[p], cones[p]);
687: maxConeSize = PetscMax(maxConeSize, cones[p]);
688: }
689: DMSetUp(dm);
690: ISGetIndices(cellsIS, &cells);
691: ISGetIndices(orntsIS, &ornts);
692: PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
693: for (p = 0, q = 0; p < pEnd; ++p) {
694: for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
695: DMPlexSetCone(dm, order[p], cone);
696: DMPlexSetConeOrientation(dm, order[p], ornt);
697: }
698: PetscFree2(cone,ornt);
699: ISRestoreIndices(orderIS, &order);
700: ISRestoreIndices(conesIS, &cones);
701: ISRestoreIndices(cellsIS, &cells);
702: ISRestoreIndices(orntsIS, &ornts);
703: ISDestroy(&orderIS);
704: ISDestroy(&conesIS);
705: ISDestroy(&cellsIS);
706: ISDestroy(&orntsIS);
707: DMPlexSymmetrize(dm);
708: DMPlexStratify(dm);
709: /* Create coordinates */
710: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
711: 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);
712: DMGetCoordinateSection(dm, &coordSection);
713: PetscSectionSetNumFields(coordSection, 1);
714: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
715: PetscSectionSetChart(coordSection, vStart, vEnd);
716: for (v = vStart; v < vEnd; ++v) {
717: PetscSectionSetDof(coordSection, v, spatialDim);
718: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
719: }
720: PetscSectionSetUp(coordSection);
721: DMSetCoordinates(dm, coordinates);
722: VecDestroy(&coordinates);
723: /* Read Labels*/
724: ctx.rank = rank;
725: ctx.dm = dm;
726: ctx.viewer = viewer;
727: PetscViewerHDF5PushGroup(viewer, "/labels");
728: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
729: PetscStackCallHDF5(H5Literate,(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx));
730: PetscStackCallHDF5(H5Gclose,(groupId));
731: PetscViewerHDF5PopGroup(viewer);
732: return(0);
733: }
734: #endif