Actual source code: plexhdf5.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc-private/isimpl.h>
3: #include <petscviewerhdf5.h>
5: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
6: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
8: #if defined(PETSC_HAVE_HDF5)
11: static PetscErrorCode GetField_Static(DM dm, PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
12: {
13: PetscInt *subIndices;
14: PetscInt Nc, subSize = 0, subOff = 0, p;
18: PetscSectionGetFieldComponents(section, field, &Nc);
19: for (p = pStart; p < pEnd; ++p) {
20: PetscInt gdof, fdof = 0;
22: PetscSectionGetDof(sectionGlobal, p, &gdof);
23: if (gdof > 0) {PetscSectionGetFieldDof(section, p, field, &fdof);}
24: subSize += fdof;
25: }
26: PetscMalloc1(subSize, &subIndices);
27: for (p = pStart; p < pEnd; ++p) {
28: PetscInt gdof, goff;
30: PetscSectionGetDof(sectionGlobal, p, &gdof);
31: if (gdof > 0) {
32: PetscInt fdof, fc, f2, poff = 0;
34: PetscSectionGetOffset(sectionGlobal, p, &goff);
35: /* Can get rid of this loop by storing field information in the global section */
36: for (f2 = 0; f2 < field; ++f2) {
37: PetscSectionGetFieldDof(section, p, f2, &fdof);
38: poff += fdof;
39: }
40: PetscSectionGetFieldDof(section, p, field, &fdof);
41: for (fc = 0; fc < fdof; ++fc, ++subOff) subIndices[subOff] = goff+poff+fc;
42: }
43: }
44: ISCreateGeneral(PetscObjectComm((PetscObject) dm), subSize, subIndices, PETSC_OWN_POINTER, is);
45: VecGetSubVector(v, *is, subv);
46: VecSetBlockSize(*subv, Nc);
47: return(0);
48: }
52: static PetscErrorCode RestoreField_Static(DM dm, PetscSection section, PetscSection sectionGlobal, Vec v, PetscInt field, PetscInt pStart, PetscInt pEnd, IS *is, Vec *subv)
53: {
57: VecRestoreSubVector(v, *is, subv);
58: ISDestroy(is);
59: return(0);
60: }
64: static PetscErrorCode DMSequenceView_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar value, PetscViewer viewer)
65: {
66: Vec stamp;
67: PetscMPIInt rank;
71: if (seqnum < 0) return(0);
72: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
73: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
74: VecSetBlockSize(stamp, 1);
75: PetscObjectSetName((PetscObject) stamp, seqname);
76: if (!rank) {
77: PetscReal timeScale;
78: PetscBool istime;
80: PetscStrncmp(seqname, "time", 5, &istime);
81: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); value *= timeScale;}
82: VecSetValue(stamp, 0, value, INSERT_VALUES);
83: }
84: VecAssemblyBegin(stamp);
85: VecAssemblyEnd(stamp);
86: PetscViewerHDF5PushGroup(viewer, "/");
87: PetscViewerHDF5SetTimestep(viewer, seqnum);
88: VecView(stamp, viewer);
89: PetscViewerHDF5PopGroup(viewer);
90: VecDestroy(&stamp);
91: return(0);
92: }
96: PetscErrorCode DMSequenceLoad_HDF5(DM dm, const char *seqname, PetscInt seqnum, PetscScalar *value, PetscViewer viewer)
97: {
98: Vec stamp;
99: PetscMPIInt rank;
103: if (seqnum < 0) return(0);
104: MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);
105: VecCreateMPI(PetscObjectComm((PetscObject) viewer), rank ? 0 : 1, 1, &stamp);
106: VecSetBlockSize(stamp, 1);
107: PetscObjectSetName((PetscObject) stamp, seqname);
108: PetscViewerHDF5PushGroup(viewer, "/");
109: PetscViewerHDF5SetTimestep(viewer, seqnum);
110: VecLoad(stamp, viewer);
111: PetscViewerHDF5PopGroup(viewer);
112: if (!rank) {
113: const PetscScalar *a;
114: PetscReal timeScale;
115: PetscBool istime;
117: VecGetArrayRead(stamp, &a);
118: *value = a[0];
119: VecRestoreArrayRead(stamp, &a);
120: PetscStrncmp(seqname, "time", 5, &istime);
121: if (istime) {DMPlexGetScale(dm, PETSC_UNIT_TIME, &timeScale); *value /= timeScale;}
122: }
123: VecDestroy(&stamp);
124: return(0);
125: }
129: PetscErrorCode VecView_Plex_Local_HDF5(Vec v, PetscViewer viewer)
130: {
131: DM dm;
132: DM dmBC;
133: PetscSection section, sectionGlobal;
134: Vec gv;
135: const char *name;
136: PetscViewerVTKFieldType ft;
137: PetscViewerFormat format;
138: PetscInt seqnum;
139: PetscReal seqval;
140: PetscBool isseq;
141: PetscErrorCode ierr;
144: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
145: VecGetDM(v, &dm);
146: DMGetDefaultSection(dm, §ion);
147: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
148: PetscViewerHDF5SetTimestep(viewer, seqnum);
149: DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer);
150: PetscViewerGetFormat(viewer, &format);
151: DMGetOutputDM(dm, &dmBC);
152: DMGetDefaultGlobalSection(dmBC, §ionGlobal);
153: DMGetGlobalVector(dmBC, &gv);
154: PetscObjectGetName((PetscObject) v, &name);
155: PetscObjectSetName((PetscObject) gv, name);
156: DMLocalToGlobalBegin(dmBC, v, INSERT_VALUES, gv);
157: DMLocalToGlobalEnd(dmBC, v, INSERT_VALUES, gv);
158: PetscObjectTypeCompare((PetscObject) gv, VECSEQ, &isseq);
159: if (isseq) {VecView_Seq(gv, viewer);}
160: else {VecView_MPI(gv, viewer);}
161: if (format == PETSC_VIEWER_HDF5_VIZ) {
162: /* Output visualization representation */
163: PetscInt numFields, f;
165: PetscSectionGetNumFields(section, &numFields);
166: for (f = 0; f < numFields; ++f) {
167: Vec subv;
168: IS is;
169: const char *fname, *fgroup;
170: char group[PETSC_MAX_PATH_LEN];
171: PetscInt pStart, pEnd;
172: PetscBool flag;
174: DMPlexGetFieldType_Internal(dm, section, f, &pStart, &pEnd, &ft);
175: fgroup = (ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_POINT_FIELD) ? "/vertex_fields" : "/cell_fields";
176: PetscSectionGetFieldName(section, f, &fname);
177: if (!fname) continue;
178: PetscViewerHDF5PushGroup(viewer, fgroup);
179: GetField_Static(dmBC, section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
180: PetscObjectSetName((PetscObject) subv, fname);
181: if (isseq) {VecView_Seq(subv, viewer);}
182: else {VecView_MPI(subv, viewer);}
183: RestoreField_Static(dmBC, section, sectionGlobal, gv, f, pStart, pEnd, &is, &subv);
184: PetscViewerHDF5PopGroup(viewer);
185: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "%s/%s", fgroup, fname);
186: PetscViewerHDF5HasAttribute(viewer, group, "vector_field_type", &flag);
187: if (!flag) {
188: if ((ft == PETSC_VTK_POINT_VECTOR_FIELD) || (ft == PETSC_VTK_CELL_VECTOR_FIELD)) {
189: PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "vector");
190: } else {
191: PetscViewerHDF5WriteAttribute(viewer, group, "vector_field_type", PETSC_STRING, "scalar");
192: }
193: }
194: }
195: }
196: DMRestoreGlobalVector(dmBC, &gv);
197: return(0);
198: }
202: PetscErrorCode VecView_Plex_HDF5(Vec v, PetscViewer viewer)
203: {
204: DM dm;
205: Vec locv;
206: const char *name;
210: VecGetDM(v, &dm);
211: DMGetLocalVector(dm, &locv);
212: PetscObjectGetName((PetscObject) v, &name);
213: PetscObjectSetName((PetscObject) locv, name);
214: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, locv);
215: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, locv);
216: PetscViewerHDF5PushGroup(viewer, "/fields");
217: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
218: VecView_Plex_Local(locv, viewer);
219: PetscViewerPopFormat(viewer);
220: PetscViewerHDF5PopGroup(viewer);
221: DMRestoreLocalVector(dm, &locv);
222: return(0);
223: }
227: PetscErrorCode VecLoad_Plex_HDF5(Vec v, PetscViewer viewer)
228: {
229: DM dm;
230: Vec locv;
231: const char *name;
232: PetscInt seqnum;
236: VecGetDM(v, &dm);
237: DMGetLocalVector(dm, &locv);
238: PetscObjectGetName((PetscObject) v, &name);
239: PetscObjectSetName((PetscObject) locv, name);
240: DMGetOutputSequenceNumber(dm, &seqnum, NULL);
241: PetscViewerHDF5SetTimestep(viewer, seqnum);
242: PetscViewerHDF5PushGroup(viewer, "/fields");
243: VecLoad_Plex_Local(locv, viewer);
244: PetscViewerHDF5PopGroup(viewer);
245: DMLocalToGlobalBegin(dm, locv, INSERT_VALUES, v);
246: DMLocalToGlobalEnd(dm, locv, INSERT_VALUES, v);
247: DMRestoreLocalVector(dm, &locv);
248: return(0);
249: }
253: static PetscErrorCode DMPlexWriteTopology_HDF5_Static(DM dm, IS globalPointNumbers, PetscViewer viewer)
254: {
255: IS orderIS, conesIS, cellsIS, orntsIS;
256: const PetscInt *gpoint;
257: PetscInt *order, *sizes, *cones, *ornts;
258: PetscInt dim, pStart, pEnd, p, conesSize = 0, cellsSize = 0, c = 0, s = 0;
259: PetscErrorCode ierr;
262: ISGetIndices(globalPointNumbers, &gpoint);
263: DMPlexGetDimension(dm, &dim);
264: DMPlexGetChart(dm, &pStart, &pEnd);
265: for (p = pStart; p < pEnd; ++p) {
266: if (gpoint[p] >= 0) {
267: PetscInt coneSize;
269: DMPlexGetConeSize(dm, p, &coneSize);
270: conesSize += 1;
271: cellsSize += coneSize;
272: }
273: }
274: PetscMalloc1(conesSize, &order);
275: PetscMalloc1(conesSize, &sizes);
276: PetscMalloc1(cellsSize, &cones);
277: PetscMalloc1(cellsSize, &ornts);
278: for (p = pStart; p < pEnd; ++p) {
279: if (gpoint[p] >= 0) {
280: const PetscInt *cone, *ornt;
281: PetscInt coneSize, cp;
283: DMPlexGetConeSize(dm, p, &coneSize);
284: DMPlexGetCone(dm, p, &cone);
285: DMPlexGetConeOrientation(dm, p, &ornt);
286: order[s] = gpoint[p];
287: sizes[s++] = coneSize;
288: for (cp = 0; cp < coneSize; ++cp, ++c) {cones[c] = gpoint[cone[cp]] < 0 ? -(gpoint[cone[cp]]+1) : gpoint[cone[cp]]; ornts[c] = ornt[cp];}
289: }
290: }
291: if (s != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of points %d != %d", s, conesSize);
292: if (c != cellsSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cone points %d != %d", c, cellsSize);
293: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, order, PETSC_OWN_POINTER, &orderIS);
294: PetscObjectSetName((PetscObject) orderIS, "order");
295: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, sizes, PETSC_OWN_POINTER, &conesIS);
296: PetscObjectSetName((PetscObject) conesIS, "cones");
297: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, cones, PETSC_OWN_POINTER, &cellsIS);
298: PetscObjectSetName((PetscObject) cellsIS, "cells");
299: ISCreateGeneral(PetscObjectComm((PetscObject) dm), cellsSize, ornts, PETSC_OWN_POINTER, &orntsIS);
300: PetscObjectSetName((PetscObject) orntsIS, "orientation");
301: PetscViewerHDF5PushGroup(viewer, "/topology");
302: ISView(orderIS, viewer);
303: ISView(conesIS, viewer);
304: ISView(cellsIS, viewer);
305: ISView(orntsIS, viewer);
306: PetscViewerHDF5PopGroup(viewer);
307: ISDestroy(&orderIS);
308: ISDestroy(&conesIS);
309: ISDestroy(&cellsIS);
310: ISDestroy(&orntsIS);
311: ISRestoreIndices(globalPointNumbers, &gpoint);
313: PetscViewerHDF5WriteAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
314: return(0);
315: }
319: static PetscErrorCode DMPlexWriteTopology_Vertices_HDF5_Static(DM dm, DMLabel label, PetscInt labelId, PetscViewer viewer)
320: {
321: IS cellIS, globalVertexNumbers;
322: const PetscInt *gvertex;
323: PetscInt *vertices;
324: PetscInt dim, depth, vStart, vEnd, v, cellHeight, cStart, cMax, cEnd, cell, conesSize = 0, numCornersLocal = 0, numCorners;
325: hid_t fileId, groupId;
326: herr_t status;
327: PetscErrorCode ierr;
330: DMPlexGetDimension(dm, &dim);
331: DMPlexGetDepth(dm, &depth);
332: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
333: DMPlexGetVTKCellHeight(dm, &cellHeight);
334: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
335: DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, PETSC_NULL);
336: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
337: if (depth == 1) return(0);
338: for (cell = cStart; cell < cEnd; ++cell) {
339: PetscInt *closure = NULL;
340: PetscInt closureSize, v, Nc = 0;
342: if (label) {
343: PetscInt value;
344: DMLabelGetValue(label, cell, &value);
345: if (value == labelId) continue;
346: }
347: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
348: for (v = 0; v < closureSize*2; v += 2) {
349: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++Nc;
350: }
351: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
352: conesSize += Nc;
353: if (!numCornersLocal) numCornersLocal = Nc;
354: else if (numCornersLocal != Nc) numCornersLocal = 1;
355: }
356: MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));
357: if (numCornersLocal && (numCornersLocal != numCorners || numCorners == 1)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Visualization topology currently only supports identical cell shapes");
359: PetscViewerHDF5PushGroup(viewer, "/viz");
360: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
361: status = H5Gclose(groupId);CHKERRQ(status);
362: PetscViewerHDF5PopGroup(viewer);
364: DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
365: ISGetIndices(globalVertexNumbers, &gvertex);
366: PetscMalloc1(conesSize, &vertices);
367: for (cell = cStart, v = 0; cell < cEnd; ++cell) {
368: PetscInt *closure = NULL;
369: PetscInt closureSize, Nc = 0, p;
371: if (label) {
372: PetscInt value;
373: DMLabelGetValue(label, cell, &value);
374: if (value == labelId) continue;
375: }
376: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
377: for (p = 0; p < closureSize*2; p += 2) {
378: if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
379: closure[Nc++] = closure[p];
380: }
381: }
382: DMPlexInvertCell_Internal(dim, Nc, closure);
383: for (p = 0; p < Nc; ++p) {
384: const PetscInt gv = gvertex[closure[p] - vStart];
385: vertices[v++] = gv < 0 ? -(gv+1) : gv;
386: }
387: DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
388: }
389: if (v != conesSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of cell vertices %d != %d", v, conesSize);
390: ISCreateGeneral(PetscObjectComm((PetscObject) dm), conesSize, vertices, PETSC_OWN_POINTER, &cellIS);
391: PetscLayoutSetBlockSize(cellIS->map, numCorners);
392: PetscObjectSetName((PetscObject) cellIS, "cells");
393: PetscViewerHDF5PushGroup(viewer, "/viz/topology");
394: ISView(cellIS, viewer);
395: #if 0
396: if (numCorners == 1) {
397: ISView(coneIS, viewer);
398: } else {
399: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
400: }
401: ISDestroy(&coneIS);
402: #else
403: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_corners", PETSC_INT, (void *) &numCorners);
404: #endif
405: ISDestroy(&cellIS);
407: PetscViewerHDF5WriteAttribute(viewer, "/viz/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
408: return(0);
409: }
413: static PetscErrorCode DMPlexWriteCoordinates_HDF5_Static(DM dm, PetscViewer viewer)
414: {
415: DM cdm;
416: Vec coordinates, newcoords;
417: PetscReal lengthScale;
421: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
422: DMGetCoordinateDM(dm, &cdm);
423: DMGetCoordinatesLocal(dm, &coordinates);
424: DMGetLocalVector(cdm, &newcoords);
425: PetscObjectSetName((PetscObject) newcoords, "vertices");
426: VecCopy(coordinates, newcoords);
427: VecScale(newcoords, lengthScale);
428: /* Use the local version to bypass the default group setting */
429: PetscViewerHDF5PushGroup(viewer, "/geometry");
430: PetscViewerPushFormat(viewer, PETSC_VIEWER_NATIVE);
431: VecView(newcoords, viewer);
432: PetscViewerPopFormat(viewer);
433: PetscViewerHDF5PopGroup(viewer);
434: DMRestoreLocalVector(cdm, &newcoords);
435: return(0);
436: }
440: static PetscErrorCode DMPlexWriteCoordinates_Vertices_HDF5_Static(DM dm, PetscViewer viewer)
441: {
442: Vec coordinates, newcoords;
443: PetscSection cSection;
444: PetscScalar *coords, *ncoords;
445: const PetscReal *L;
446: PetscReal lengthScale;
447: PetscInt vStart, vEnd, v, bs, coordSize, dof, off, d;
448: hid_t fileId, groupId;
449: herr_t status;
450: PetscErrorCode ierr;
453: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
454: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
455: DMGetCoordinatesLocal(dm, &coordinates);
456: VecGetBlockSize(coordinates, &bs);
457: VecGetLocalSize(coordinates, &coordSize);
458: if (coordSize == (vEnd - vStart)*bs) return(0);
459: DMGetPeriodicity(dm, NULL, &L);
460: DMGetCoordinateSection(dm, &cSection);
461: VecCreate(PetscObjectComm((PetscObject) coordinates), &newcoords);
462: coordSize = 0;
463: for (v = vStart; v < vEnd; ++v) {
464: PetscSectionGetDof(cSection, v, &dof);
465: if (L && dof == 2) coordSize += dof+1;
466: else coordSize += dof;
467: }
468: if (L && bs == 2) {VecSetBlockSize(newcoords, bs+1);}
469: else {VecSetBlockSize(newcoords, bs);}
470: VecSetSizes(newcoords, coordSize, PETSC_DETERMINE);
471: VecSetType(newcoords,VECSTANDARD);
472: VecGetArray(coordinates, &coords);
473: VecGetArray(newcoords, &ncoords);
474: coordSize = 0;
475: for (v = vStart; v < vEnd; ++v) {
476: PetscSectionGetDof(cSection, v, &dof);
477: PetscSectionGetOffset(cSection, v, &off);
478: if (L && dof == 2) {
479: ncoords[coordSize++] = -cos(2.0*PETSC_PI*coords[off+0]/L[0])*(L[0]/(2.0*PETSC_PI));
480: ncoords[coordSize++] = coords[off+1];
481: ncoords[coordSize++] = sin(2.0*PETSC_PI*coords[off+0]/L[0])*(L[0]/(2.0*PETSC_PI));
482: } else {
483: for (d = 0; d < dof; ++d, ++coordSize) ncoords[coordSize] = coords[off+d];
484: }
485: }
486: VecRestoreArray(coordinates, &coords);
487: VecRestoreArray(newcoords, &ncoords);
488: PetscObjectSetName((PetscObject) newcoords, "vertices");
489: VecScale(newcoords, lengthScale);
490: /* Use the local version to bypass the default group setting */
491: PetscViewerHDF5PushGroup(viewer, "/viz");
492: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
493: status = H5Gclose(groupId);CHKERRQ(status);
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: herr_t status;
514: PetscViewerFormat format;
515: PetscErrorCode ierr;
518: DMPlexCreatePointNumbering(dm, &globalPointNumbers);
519: PetscViewerGetFormat(viewer, &format);
520: DMPlexWriteCoordinates_HDF5_Static(dm, viewer);
521: if (format == PETSC_VIEWER_HDF5_VIZ) {DMPlexWriteCoordinates_Vertices_HDF5_Static(dm, viewer);}
522: DMPlexWriteTopology_HDF5_Static(dm, globalPointNumbers, viewer);
523: if (format == PETSC_VIEWER_HDF5_VIZ) {DMPlexWriteTopology_Vertices_HDF5_Static(dm, label, labelId, viewer);}
524: /* Write Labels*/
525: ISGetIndices(globalPointNumbers, &gpoint);
526: PetscViewerHDF5PushGroup(viewer, "/labels");
527: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
528: if (groupId != fileId) {status = H5Gclose(groupId);CHKERRQ(status);}
529: PetscViewerHDF5PopGroup(viewer);
530: DMPlexGetNumLabels(dm, &numLabels);
531: for (l = 0; l < numLabels; ++l) {
532: DMLabel label;
533: const char *name;
534: IS valueIS, globalValueIS;
535: const PetscInt *values;
536: PetscInt numValues, v;
537: PetscBool isDepth;
538: char group[PETSC_MAX_PATH_LEN];
540: DMPlexGetLabelName(dm, l, &name);
541: PetscStrncmp(name, "depth", 10, &isDepth);
542: if (isDepth) continue;
543: DMPlexGetLabel(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) {status = H5Gclose(groupId);CHKERRQ(status);}
548: PetscViewerHDF5PopGroup(viewer);
549: DMLabelGetValueIS(label, &valueIS);
550: ISAllGather(valueIS, &globalValueIS);
551: ISSortRemoveDups(globalValueIS);
552: ISGetLocalSize(globalValueIS, &numValues);
553: ISGetIndices(globalValueIS, &values);
554: for (v = 0; v < numValues; ++v) {
555: IS stratumIS, globalStratumIS;
556: const PetscInt *spoints;
557: PetscInt *gspoints, n, gn, p;
558: const char *iname;
560: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%d", name, values[v]);
561: DMLabelGetStratumIS(label, values[v], &stratumIS);
563: ISGetLocalSize(stratumIS, &n);
564: ISGetIndices(stratumIS, &spoints);
565: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) ++gn;
566: PetscMalloc1(gn,&gspoints);
567: for (gn = 0, p = 0; p < n; ++p) if (gpoint[spoints[p]] >= 0) gspoints[gn++] = gpoint[spoints[p]];
568: ISRestoreIndices(stratumIS, &spoints);
569: ISCreateGeneral(PetscObjectComm((PetscObject) dm), gn, gspoints, PETSC_OWN_POINTER, &globalStratumIS);
570: PetscObjectGetName((PetscObject) stratumIS, &iname);
571: PetscObjectSetName((PetscObject) globalStratumIS, iname);
573: PetscViewerHDF5PushGroup(viewer, group);
574: ISView(globalStratumIS, viewer);
575: PetscViewerHDF5PopGroup(viewer);
576: ISDestroy(&globalStratumIS);
577: ISDestroy(&stratumIS);
578: }
579: ISRestoreIndices(globalValueIS, &values);
580: ISDestroy(&globalValueIS);
581: ISDestroy(&valueIS);
582: }
583: ISRestoreIndices(globalPointNumbers, &gpoint);
584: ISDestroy(&globalPointNumbers);
585: return(0);
586: }
588: typedef struct {
589: PetscMPIInt rank;
590: DM dm;
591: PetscViewer viewer;
592: DMLabel label;
593: } LabelCtx;
595: static herr_t ReadLabelStratumHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
596: {
597: PetscViewer viewer = ((LabelCtx *) op_data)->viewer;
598: DMLabel label = ((LabelCtx *) op_data)->label;
599: IS stratumIS;
600: const PetscInt *ind;
601: PetscInt value, N, i;
602: const char *lname;
603: char group[PETSC_MAX_PATH_LEN];
604: PetscErrorCode ierr;
606: PetscOptionsStringToInt(name, &value);
607: ISCreate(PetscObjectComm((PetscObject) viewer), &stratumIS);
608: PetscObjectSetName((PetscObject) stratumIS, "indices");
609: DMLabelGetName(label, &lname);
610: PetscSNPrintf(group, PETSC_MAX_PATH_LEN, "/labels/%s/%s", lname, name);
611: PetscViewerHDF5PushGroup(viewer, group);
612: {
613: /* Force serial load */
614: PetscViewerHDF5ReadSizes(viewer, "indices", NULL, &N);
615: PetscLayoutSetLocalSize(stratumIS->map, !((LabelCtx *) op_data)->rank ? N : 0);
616: PetscLayoutSetSize(stratumIS->map, N);
617: }
618: ISLoad(stratumIS, viewer);
619: PetscViewerHDF5PopGroup(viewer);
620: ISGetLocalSize(stratumIS, &N);
621: ISGetIndices(stratumIS, &ind);
622: for (i = 0; i < N; ++i) {DMLabelSetValue(label, ind[i], value);}
623: ISRestoreIndices(stratumIS, &ind);
624: ISDestroy(&stratumIS);
625: return 0;
626: }
628: static herr_t ReadLabelHDF5_Static(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
629: {
630: DM dm = ((LabelCtx *) op_data)->dm;
631: hsize_t idx;
632: herr_t status;
635: DMPlexCreateLabel(dm, name); if (ierr) return (herr_t) ierr;
636: DMPlexGetLabel(dm, name, &((LabelCtx *) op_data)->label); if (ierr) return (herr_t) ierr;
637: status = H5Literate_by_name(g_id, name, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelStratumHDF5_Static, op_data, 0);
638: return status;
639: }
643: /* The first version will read everything onto proc 0, letting the user distribute
644: The next will create a naive partition, and then rebalance after reading
645: */
646: PetscErrorCode DMPlexLoad_HDF5(DM dm, PetscViewer viewer)
647: {
648: LabelCtx ctx;
649: PetscSection coordSection;
650: Vec coordinates;
651: IS orderIS, conesIS, cellsIS, orntsIS;
652: const PetscInt *order, *cones, *cells, *ornts;
653: PetscReal lengthScale;
654: PetscInt *cone, *ornt;
655: PetscInt dim, spatialDim, N, numVertices, vStart, vEnd, v, pEnd, p, q, maxConeSize = 0, c;
656: hid_t fileId, groupId;
657: hsize_t idx = 0;
658: herr_t status;
659: PetscMPIInt rank;
660: PetscErrorCode ierr;
663: MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
664: /* Read toplogy */
665: PetscViewerHDF5ReadAttribute(viewer, "/topology/cells", "cell_dim", PETSC_INT, (void *) &dim);
666: DMPlexSetDimension(dm, dim);
667: PetscViewerHDF5PushGroup(viewer, "/topology");
669: ISCreate(PetscObjectComm((PetscObject) dm), &orderIS);
670: PetscObjectSetName((PetscObject) orderIS, "order");
671: ISCreate(PetscObjectComm((PetscObject) dm), &conesIS);
672: PetscObjectSetName((PetscObject) conesIS, "cones");
673: ISCreate(PetscObjectComm((PetscObject) dm), &cellsIS);
674: PetscObjectSetName((PetscObject) cellsIS, "cells");
675: ISCreate(PetscObjectComm((PetscObject) dm), &orntsIS);
676: PetscObjectSetName((PetscObject) orntsIS, "orientation");
677: {
678: /* Force serial load */
679: PetscViewerHDF5ReadSizes(viewer, "order", NULL, &pEnd);
680: PetscLayoutSetLocalSize(orderIS->map, !rank ? pEnd : 0);
681: PetscLayoutSetSize(orderIS->map, pEnd);
682: PetscViewerHDF5ReadSizes(viewer, "cones", NULL, &pEnd);
683: PetscLayoutSetLocalSize(conesIS->map, !rank ? pEnd : 0);
684: PetscLayoutSetSize(conesIS->map, pEnd);
685: PetscViewerHDF5ReadSizes(viewer, "cells", NULL, &N);
686: PetscLayoutSetLocalSize(cellsIS->map, !rank ? N : 0);
687: PetscLayoutSetSize(cellsIS->map, N);
688: PetscViewerHDF5ReadSizes(viewer, "orientation", NULL, &N);
689: PetscLayoutSetLocalSize(orntsIS->map, !rank ? N : 0);
690: PetscLayoutSetSize(orntsIS->map, N);
691: }
692: ISLoad(orderIS, viewer);
693: ISLoad(conesIS, viewer);
694: ISLoad(cellsIS, viewer);
695: ISLoad(orntsIS, viewer);
696: PetscViewerHDF5PopGroup(viewer);
697: /* Read geometry */
698: PetscViewerHDF5PushGroup(viewer, "/geometry");
699: VecCreate(PetscObjectComm((PetscObject) dm), &coordinates);
700: PetscObjectSetName((PetscObject) coordinates, "vertices");
701: {
702: /* Force serial load */
703: PetscViewerHDF5ReadSizes(viewer, "vertices", &spatialDim, &N);
704: VecSetSizes(coordinates, !rank ? N : 0, N);
705: VecSetBlockSize(coordinates, spatialDim);
706: }
707: VecLoad(coordinates, viewer);
708: PetscViewerHDF5PopGroup(viewer);
709: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
710: VecScale(coordinates, 1.0/lengthScale);
711: VecGetLocalSize(coordinates, &numVertices);
712: VecGetBlockSize(coordinates, &spatialDim);
713: numVertices /= spatialDim;
714: /* Create Plex */
715: DMPlexSetChart(dm, 0, pEnd);
716: ISGetIndices(orderIS, &order);
717: ISGetIndices(conesIS, &cones);
718: for (p = 0; p < pEnd; ++p) {
719: DMPlexSetConeSize(dm, order[p], cones[p]);
720: maxConeSize = PetscMax(maxConeSize, cones[p]);
721: }
722: DMSetUp(dm);
723: ISGetIndices(cellsIS, &cells);
724: ISGetIndices(orntsIS, &ornts);
725: PetscMalloc2(maxConeSize,&cone,maxConeSize,&ornt);
726: for (p = 0, q = 0; p < pEnd; ++p) {
727: for (c = 0; c < cones[p]; ++c, ++q) {cone[c] = cells[q]; ornt[c] = ornts[q];}
728: DMPlexSetCone(dm, order[p], cone);
729: DMPlexSetConeOrientation(dm, order[p], ornt);
730: }
731: PetscFree2(cone,ornt);
732: ISRestoreIndices(orderIS, &order);
733: ISRestoreIndices(conesIS, &cones);
734: ISRestoreIndices(cellsIS, &cells);
735: ISRestoreIndices(orntsIS, &ornts);
736: ISDestroy(&orderIS);
737: ISDestroy(&conesIS);
738: ISDestroy(&cellsIS);
739: ISDestroy(&orntsIS);
740: DMPlexSymmetrize(dm);
741: DMPlexStratify(dm);
742: /* Create coordinates */
743: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
744: 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);
745: DMGetCoordinateSection(dm, &coordSection);
746: PetscSectionSetNumFields(coordSection, 1);
747: PetscSectionSetFieldComponents(coordSection, 0, spatialDim);
748: PetscSectionSetChart(coordSection, vStart, vEnd);
749: for (v = vStart; v < vEnd; ++v) {
750: PetscSectionSetDof(coordSection, v, spatialDim);
751: PetscSectionSetFieldDof(coordSection, v, 0, spatialDim);
752: }
753: PetscSectionSetUp(coordSection);
754: DMSetCoordinates(dm, coordinates);
755: VecDestroy(&coordinates);
756: /* Read Labels*/
757: ctx.rank = rank;
758: ctx.dm = dm;
759: ctx.viewer = viewer;
760: PetscViewerHDF5PushGroup(viewer, "/labels");
761: PetscViewerHDF5OpenGroup(viewer, &fileId, &groupId);
762: status = H5Literate(groupId, H5_INDEX_NAME, H5_ITER_NATIVE, &idx, ReadLabelHDF5_Static, &ctx);CHKERRQ(status);
763: status = H5Gclose(groupId);CHKERRQ(status);
764: PetscViewerHDF5PopGroup(viewer);
765: return(0);
766: }
767: #endif